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Determination of optimum hull form 
for passenger car ferry with regard 
to its sea-keeping qualities and additional 
resistance in waves 


Tomasz Cepowski, Ph. D. 
Szczecin Maritime University 


ABSTRACT 


This paper presents a method which makes it possible to determine optimum hull form of 
passenger car ferry with regard to selected sea-keeping qualities and additional resistance 
in waves. In the first phase of investigations a hull form characterized by the highest 


y qualities was selected from the list of similar ships. Next, its optimum dimension ratios 


t 


were determined. Design criteria were formulated by using a method based on deterministic 
‘ scenarios, but objective functions of partial targets were determined in the form of artificial 
neural networks. To select the best design variants elements of fuzzy logic were used, that 


made it possible a.o. to show merits of the design by means of linguistic variables. Such approach made 
it possible to find the best hull form and its dimensions from the point of view of all considered criteria 
simultaneously. 


Keywords: sea-keeping qualities, roll-on-roll-off ferry ship, rolling, transverse accelerations, vertical accelerations, motion 
sickness index, additional resistance in waves, ship hull form, ship design parameters, artificial neural networks, fuzzy logic. 


INTRODUCTION 


In the ship design process design solutions complying both 
with economic criteria and technical constraints, are searched 
for. The economic criteria contain a set of ship owner’s 
requirements to which can belong a.o. ship service speed of 
a significant impact on ship operation profitability on a given 
shipping route. As the ship is often operated in heavy weather 
conditions the reaching of its assumed service speed depends 
a.o. on its additional resistance in waves. And, for certain types 
of ships, e.g. passenger car ferries , their non-susceptibility 
to weather conditions, i.e. the so called good sea-keeping 
qualities, constitutes a more - and- more important technical 
limitation. The additional ship resistance in waves and good 
sea-keeping qualities are strongly dependent on ship’s hull 
form and dimensions. Therefore modeling the ship qualities 
should be performed already in the preliminary ship design 
stage. Selection of a non-suitable hull form and dimensions of 
a ship irreversibly worsens its design merits and any change 
of dimensions of a built ship is economically unprofitable. In 
[1, 3] a method which makes it possible to take into account 
the ship sea-keeping qualities and additional ship resistance 
in waves in the phase of parametric design, is presented. As 
a result of the research, optimum values of ship dimensions and 
certain global coefficients characterizing hull form of a ro-ro 
ferry ship were determined. However the aim of the research 
presented in this paper was to model a geometrically described 
hull form and to determine an optimum hull form with regard 


to its sea-keeping qualities and additional resistance in waves. 
The obtained results are presented with the use of an example 
design task. 


DESIGN TASK 


In the research the following design tasks were formulated: 
to find such hull form of the ro-ro ferry, which minimizes values 
of the design criteria at fulfilled assumed constraints for the 
passenger car ferry sailing in heavy weather conditions. 

In the design task the following were assumed: 


\? 


** ship design parameters describing its hull form, considered 
independent variables: 


L/B 


L, B, d—ship length, breadth and draught, respectively, 
CB — block coefficient of underwater part of hull, 
CWL - waterplane coefficient; 


% constraints: 

theoretical displacement V = 13723 m° 

range of L/B : from 5.17 to 6.74 

range of B/d : from 3.22 to 4.46 

range of CB : from 0.6 to 0.64 

range of CWL: from 0.8 to 0.85 

the initial transverse metacentric height GM = 1 m 


+ 


++ +% + o 
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METHOD 


To solve the problem was elaborated a method consisted 
of the two phases: 
1. On the basis of the list of similar ships - determination of 
alternative design variants and selection of the best one 
(Fig. 1) 


Design assumptions 


V = const 


List of alternative variant 
® Cg =var 
e L/B, B/d = var 


Deterministic 
scenarios 


Calculation of values 
of design criteria by means 
of exact methods 


Design 
criteria 


Choice of the best variant 
of determined form 
CB = const = a 
CWL = const = B 


Fig. 1. Algorithm showing the 1" phase of the research 


2. On the basis of the best hull form variant selected in the 
1* phase — selection of the optimum ratios of ferry ship 
dimensions (Fig. 2). 


List of model design variants 
e CB=a 
e CWL=8 Deterministic 
e LB. Bid = var scenarios 
Calculation of model : 
values of design criteria (= Design 
by using exact methods ee 
Approximation of model 
values of design criteria 
with regard to L/B and B/d 
Determination of optimum 
values of L/B and B/d Assessment 
with regard to the best a by using 
design merits fuzzy logic 


dh 


Verification of design merits 
of the optimum variant 
by using exact 
numerical methods 


Fig. 2. Algorithm showing the 2" phase of the research 
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The assessment process of design variants carried out in 
both research phases was based on the below given design 
criteria. 


DESIGN CRITERIA 


One of the crucial design problems of ro-ro ferries is to 
ensure appropriate damage stability. The ship stability criteria 
being in force are rather far from taking into account the complex 
physical phenomena associated with ship hydromechanics, and 
they are based first of all on physical parameters calculated for 
a ship in still water conditions. Hence the application of design 
criteria inadequate to the conditions in which the ship operates 
does not ensure a sufficient level of safety. The IMO has 
initiated efforts aimed at undertaking a research on elaboration 
of new standards for ship stability assessment. One of the ways 
to solve the problem is elaboration of methods for investigating 
the intensity of ship rolling motion (e.g. [12]) or probability 
of ship capsizing (e.g. [10]). In the case of ro-ro ferry ships 
such investigations are aimed at elaboration and improvement 
of simulation models of sinking process of damaged ferry in 
given conditions (acc. [8]). The example models are presented 
in [13, 15]. The above referred research makes it possible to 
simulate ferry sinking process however any design criteria have 
not yet evolved from that. 

The design criteria assumed in this research are to an 
extent associated with the above mentioned reccomendations, 
as thy are based on statistical values of ship response in 
waves. For their formulation a method based on deterministic 
scenarios [16], was used, hence the following scenarios were 
assumed: 


Scenario no. 1: the design ferry sails with the speed v = 15 
kn in the statistical heavy weather conditions (the significant 
wave height of 3 m, the characteristic wave period in the range 
from 3s to 18 s.), the ship encounters the wave from the most 
unfavourable directions. 


Scenario no. 2: the designed ferry is in the following 
damage state : its side plating amidships has been broken at 
the level of car deck, the disable ship is in the statistical heavy 
weather conditions (the significant wave height of 3 m, the 
characteristic wave period in the range from 3s to 18 s.), the 
ship heads the wave at the encounter angle B = 120° (where: 
180° — head wave, 0° — following wave). 

Such approach made it possible to formulate the following 
design criteria: 

For Scenario no. 1: 

1. As- small- as- possible value of the motion sickness index 
MSI (acc. ISO 2631/3), [9] for the wave encouter angle B 
= 120° 

2. As- small- as- possible value of the additional ship resistance 
in waves R, for B = 180° 

3. As-—small-as-possible significant value of the roll angle @ 
, for B = 30° 

4. As-small-as-possible value of the transverse accelerations 
at the car deck, a, , acc. [16], in the point P2 in Fig. 3 , for 
B = 30° 

5. As-small-as-possible value of the vertical accelerations at 
the car deck, a, , acc. [16], in the point P2 in Fig. 3. 


For Scenario no. 2: 
6. As-low-as-possible number of events of flooding the car 
deck in the place of broken side plating, L,,,,, for the wave 


ancounter angle B = 120°, in the point P1 in Fig. 3. 


Fulfilment of the criteria no. 1, 3, 5 ensures an appropriate 
voyage comfort for passengers , moreover the criterion no. 3 


ensures safety of ship stability, resulting from a lower intesity 
of ship rolling motions. Fulfilment of the criterion no. 4 leads 
to a lower hazard due to possible shift of vehicles on the car 
deck. Fulfilment of the criterion no. 6 can lower a number of 
flooding events of the car deck due to broken side plating, and 
thus decrease mass of the water flooded on the deck , that can 
make sinking time of the ship longer. And, fulfilment of the 
criterion no. 2 can lower a decrease of ship speed in waves, and 
thus operational costs due to lower fuel consumption. 


Fig. 3. The coordinates of the points P1, P2, P3 for which the following 
quantities were calculated : number of car deck flooding events — in the 
point P1 (0.5L, 0.49B, d+Im), transverse accelerations — in the point P2 
(0.75L, 0.45B, d), vertical accelerations — in the point P3 (0.8L, 0.32B, 2d), 
where: L — ship length, B — ship breadth, d — ship draught 


THE 1* PHASE OF THE RESEARCH 


The 1* phase of the research was aimed at the finding - on 
the basis of the list of alternative variants — a hull form of the 
best merits. In the investigations were used eight ro-ro ferries 
whose dimensions and hull forms are presented in Tab. | and 
Fig. 4. 


Tab. 1. List of alternative design variants, V = 13723 m’, CB, L/B, B/d = var 


Design variants 
Variant no. 
Variant no. 
Variant no. 


Variant no. 
Variant no. 
Variant no. 
Variant no. 
Variant no. 


For each of the design variants given in Tab. | statistical 
values of the assumed design criteria were calculated. 
The conventional ITTC wave spectrum was assumed. The 
calculations were carried out by means of the exact numerical 
methods included in the SEAWAY software based on the theory 
of two-dimensional flow. For calculation of hydrodynamic 
coefficients the Frank’s method was applied [5]. The accuracy 
tests of the SEAWAY software, presented in [9], have shown 
that its calculation accuracy is high. For calculations of the 
additional ship resistance in waves the Gerritsma-Beukelman 
method was used [6]. 

To the investigations were applied such values of the design 
criteria which reached their maxima in the assumed statistical 
wave conditions (i.e. the significant wave height equal to 3 m, 
the characteristic wave period within the range from 3s to 18 s.). 
It enabled to eleminate influence of the wave parameters on ship 
responses and to lower the number of independent variables 
in the design task in question. The results of calculations are 
presented in Tab. 2. 


Variant no. 1 


nmr 
AE [> 
UUL 


E R 


EM 


BE yM 


Variant no. 


apih dPR 
7 j 


E *, $3 
CE st / Pt 
S ae $ 97 d/p 
s ioe 5 arena nae naga 


Variant no. 6 


sa 
Son 


Variant no. 7 


Variant no. 8 
Fig. 4. Hull forms of alternative design variants of ro-ro ferry, V = 13723 m, CB, L/B, B/d = var 
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Tab. 2. Values of the design criteria used in the 1” phase of the research, 
where: MSI — motion sickness index, a,— transverse acceleration at the car 
deck, a,— vertical accelerations, o — significant amplitude of ship rolling, 
R — additional ship resistance in waves, L,,,,— hourly number of car deck 
flooding events in the place of broken side plating 


Design 
variants 


Variant no. 


Variant no. 


Variant no. 


Variant no. 


Variant no. 


Variant no. 


Variant no. 


Variant no. 


In the further part of the research the best design variant 
was selected with the use of a method based on fuzzy logic 
elements. 


APPLICATION OF FUZZY LOGIC 
ELEMENTS TO SELECTION 
OF THE BEST DESIGN SOLUTION 


Fuzzy logic finds its application wherever the use of classic 
logic creates a problem consisted in difficulty of mathematical 
description of a process or when to calculate or select variables 
for solving a problem is not possible. In the design problem in 
question fuzzy logic was used as an element which aids to assess 
quality of a design variant (of low, mean or high merits) and to 
decide on choice of the best design (of the best merits). 

In the 1“ research phase the set of solutions given in Tab. 2 
was fuzzyfied with the use of the attribution function decribed 
by the Eq. (1): 


U(x,) = 0 for Xi Sig) (1) 
Xi — Ay) 


Qo) T Ay) 


where: 

u(x)  — value of attribution function for i-th design criterion 
(i= MSI, 9, a, a, R, Low) 

X — value of i-th design criterion (from Tab. 2) 


i 


Aia ba — threshold values acc. Tab. 3. 


The graphical form of the function (1) as well as its 
liguistic interpretation is presented in Fig. 5; and, values of the 


coefficients a; O are given in Tab. 3. 
a 20) 


Tab. 3. Values of the coefficients a, and a, for the design criteria: 
MSI — motion sickness index, a , ~ transverse acceleration at the car deck, 
a,— vertical accelerations, » — significant amplitude of ship rolling, 


R — additional ship resistance in waves, L,,,,— hourly number 


of car deck flooding events in the sector of broken side plating 


MSI a, a, R Lew 
[%] [m/s?] | [m/s*] | [KN] | [1/h] 


164.13 | 274.97 
179.97 | 282.43 
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LCX;) 


Design of high merits 


Design of mean merits 


Design of low merits 


OG) l 


Aia) 
Fig. 5. Attribution function acc. Eq. (1) 


Next, assessment of the design solutions fuzzyfied with 
the use of Eq. (1) was performed. To assess design quality 
the minimizing operator MIN was applied in accordance with 


Eq. (2): 
MIN = MIN(u(MSD), u(R), u(ọ), u(a), u(a ), WL ew) 
(2) 


where: 

MIN — minimizing operator 

u(MSI), w(R), u(ọ), u(a), (a ), W(L,,,) — attribution function 
of the design criteria: MSI — motion sickness index, 
a, — transverse acceleration at the car deck, a, — vertical 
accelerations, ọ — significant amplitude of ship rolling, 
R — additional ship resistance in waves, Lįẹ„ — hourly 
number of car deck flooding events in the place of broken 
side plating, respectively. 


The values of the MIN operator as well as attribution 
function of the design criteria are given in Tab. 4 and Fig. 6. 


Tab. 4. The values of the MIN operator as well as attribution function of the 
design criteria: MSI — motion sickness index, a , ~ transverse acceleration 
at the car deck, a, — vertical accelerations, o — significant amplitude of ship 
rolling, R — additional ship resistance in waves, L,— hourly number of car 
deck flooding events in the place of broken side plating 


Variants |u(MST) 


From Tab. 4 and Fig. 6 it results that variant no. 8 is of the 
best merits. 


Fig. 6. Values of the MIN operator and attribution functions 
of the design variants 


|] apeusn 
ira au 
Peay o u(av 
l lE ae 
0.9 E] m (R) 
0.8 E mu(LGW) 
0.7 ei EMIN 
0.6 | 


~ 
E 
3 


Variant 6 
Variant 5 
Variant 4 
Variant 3 
Variant 2 
Variant 1 


Fig. 6. Values of the MIN operator and attribution functions 
of the design variants 


THE 2" PHASE OF THE RESEARCH 


The 2™ phase was aimed at determination of optimum 
design parameters of the ro-ro ferry under the following 


assumptions: 
— L/B=var 
— B/d=var 


at the determined hull form complying with Variant no. 8, 
kes 

— CB =0.64 

— CWL = 0.801 


as well as the earlier accepted design assumptions. 


In the 1* step a model set of design variants was elaborated 
and values of the design criteria were calculated. The same 
assumptions, calculation method and design criteria as in the 
1* phase, were used. The calculation results are presented in 
Tab. 5. 

Next, approximations of values of the dependent variables 
were elaborated depending on the independent variables. The 


Tab. 5. The model set of independent variables and dependent ones assumed in the design task of CB = 0,64, CWL = 0,801, where: L — ship length, B — ship 
breadth, d — ship draught, MSI — motion sickness index, a,— transverse acceleration at the car deck, a, — vertical accelerations, ~ — significant amplitude of 
ship rolling, R — additional ship resistance in waves, L_,,,— hourly number of car deck flooding events in the place of broken side plating 


Independent 
variables 


Dependent variables 


Variants L/B Bla 


I- Bie ee 


Variant no. 1 


a, [m/s?] Low [1h] 


a, 
[m/s?] 


Variant no.2 


Variant no.3 


Variant no.4 


Variant no.5 


Variant no.6 


Variant no.7 


Variant no.8 


Variant no.9 


Variant no.10 


Variant no. 11 


Variant no.12 


Variant no.13 


Variant no.14 


Variant no.15 


Variant no.16 


Variant no.17 


Variant no.18 


Variant no.19 


Variant no.20 


Variant no.21 


Variant no.22 


Variant no.23 


Variant no.24 


Variant no.25 
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approximations were obtained by means of the artificial neural networks of MLP type (Multi Layer Perceptron) and presented 
with the use of Eqs. (3), (4), (5), (6), (7) and (8): 


1 
(UB B/dxS+PKA—B) X CT o |7 OL 6) 


a 
Qh 


where: 
MSI - motion sickness index [%]; L— ship length; B — ship breadth; d — ship draught; 
A — matrix of weighing factors: 


[ -0.311 0.921 -0.694 5.134 2.302 -1.432 -3.048 -3.136 0.695 0.098 -0.090 ] 
| 0.957 2.041 0.917 2.913 1.206 -0.194 -0.191 -1.058 0.459 -0.254 0.123 | 


S — matrix of coefficients: 


0.284 0 
L 0 0.806 _| 


B — vector of threshold values: 
[-0.271 -1.606 -0.124 -1.718 -0.701 0.526 0.322 -3.638 -0.342 -0.260 -0.704 ] 
C — column vector of weighing factors: 
[ -0.955 -1.137 -0.771 2.374 0.855 -1.171 -1.886 1.953 0.008 -0.537 -0.535 ] 
P — vector of translation values: [ -0.918 -2.597 | 
— coefficients of the values: a, = 0.8, a, = -0.245, a, = 0.031 


1 
J4 eB .BidxS+P)xA-B) 
R= 


QO), A, O, 


xC |e (4) 


Q2 


where: 
R — additional ship resistance in waves, [kN] 
A — matrix of values of weighing factors: 


-0.533 -2.682 4.015 0.886 -3.133 | 
0.065 -0.252 1.285 1.316 -0.798 | 
B — vector of threshold values: [-0.466 0.249 2.501 1.186 -2.844] 


C — column vector of values of weighing factors: [-0.350 -2.402 -1.706 1.705 1.517] 
ay 0,, a, — coefficient of the values: a, = -0.32, a, = -1.846, a, = 0.013 


1 
14+ e{IL/B.B/d]xS+P)xA-B) xC-a | - s 
ọ= 


0) 


where: 
ọ — significant roll amplitude [°] 
A — matrix of values of weighing factors: 


3.869 5.902 
0.798 0.478 
B — vector of threshold values: [ 0.530 2.873 ] 


C — column vector of values of weighing factors: [3.143 -2.295 ] 
ay @,, a, — coefficients of the values: a, = 0.798, a, = -0.677, a, = 0.123 


1 
xC-a,]-a 
q+ Bid kS+P)xA-B) .) 1 (6) 


a = 
v a 
2 


where: 
a, — vertical accelerations [m/s?], 
A — matrix of values of weighing factors: 


| 1.76 038 2.19 0.58 -0.92 0.89 0.19 0.50 -0.21 -0.14 0.70 0.88 0.11 0.47 3.88 -1.57 | 
| 0.12 -0.40 0.12 1.02 0.25 -0.84 0.98 -0.29 -1.12 -0.14 0.95 -0.30 0.48 -0.11 1.18 -0.32 | 
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B — vector of threshold values: 


[-0.85 -0.71 -0.62 0.15 0.86 -0.27 -0.61 


-0.89 0.34 0.05 0.60 0.17 0.39 -0.32 0.60 -1.51] 


C — column vector of values of weighing factors: 


[0.78 -0.48 1.02 -0.46 -0.77 -0.61 


-1.54 -0.41 0.79 0.08 


-1.62 -0.73 -0.46 -0.51 2.74 1.75] 


Qp Q, a, — coefficients of the values: a, = 0.157, a, = -3.727, a, = 2.272 


1 
J4 e MLB Bld xS+P)xA-B) xC a) -a (7) 
t= CA 
where: 
a, — transverse accelerations at the car deck, [m/s?] 
A — matrix of values of weighing factors: 
| 3.862 -6.941 | 
| 1.672 -1.839 | 
B — vector of threshold values: [ 2.487 -4.718 ] 
C — column vector of values of weighing factors: [2.662 2.562 | 
ay O,, a, — coefficients of the values : a, = 2.449, a, = -0.674, a, = 0.694 
1 
J4 eo WLB.BidxS+P)xA-B) xC a) OQ; (8) 
L Z +e 
GW Gs 
where: 


a, — vertical accelerations, [m/s”] 
A — matrix of values of weighing factors: 


| 0.108 -0.507 -0.278 -0.658 2.414 
| 0.363 0.154 0.064 0.538 0.079 


-3.652 


-1.223 
-0.409 -0.370 0.887 -0.521 


-0.748 -0.307 1.202 0.724 0.666 | 
-0.856 0.533 0.346 -0.174 0.424 | 


1.756 0.941 


B — vector of threshold values: 


[ -1.007 0.932 0.823 0.421 -1.062 


-0.243 0.616 -0.895 0.089 -0.700 -0.185 0.486 -0.241 


-0.286 ] 


C — column vector of values of weighing factors: 


[ -0.781 0.436 0.698 0.581 1.399 


-2.512 -0.344 1.135 -0.733 0.997 0.285 -1.148 -0.661 


-0.698 ] 


Q,, Q, a, — coefficients of the values: a, = -0.446, a, = -4.874, a, = 0.021 


In Tab. 6. the above mentioned networks are described and 
their selected statistical parameters are given. From the data 
it results that the functions (3), (4), (5), (6), (7) and (8) have 
a simple structure and high accuracy. 


Tab. 6. Structure and selected statistical parameters 
of the elaborated neural networks 


Mean Square error 


Neural RMS 


network 
structure 


Correlation 


Variable 


coefficient R ; of 
of learning 


validation 
2x11x1 
2x2x1 
2x5x1 
2x2x1 
2x16x1 
2x14x1 


In the further part of the research a hull form of the 
passenger car ferry was searched for with regard to the assumed 
criteria. To solve the problem the functions (2), (3), (4), (5), 
(6) and (7) were used. 


In searching for a design solution which satisfies the 
assumed criteria the methods of one —criterial and multi-criterial 
optimization can be applied. However the main drawback of the 
methods is that they inform only where an optimum is located 
and they do not provide any information on a form of objective 
function, hence there is no certainty that the obtained solution 
is global optimum [4]. Number of possible combinations is 
relatively small as in the design task in question only two 
dependent variables of a rather small variability range appear. 
Therefore in order to determine an optimum solution the set 
of all possible design solutions was taken into account. For 
calculations the following increments of the design parameters 
were assumed: 
> AL/B=0.01 
> AB/d=0.01 


The calculation results were analyzed in the same way as 
in the 1“ phase of the research. To this end all design solutions 
were fuzzyfied by means of Eq. (1) and for their assessment the 
minimizing operator MIN acc. Eq. (2) was applied. 


From the performed investigations it results that the best 
solution was obtaind for the following values of the design 
parameters : 
> L/B=6.72 + 6.75 
> B/d=3.57 +3.80 
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The value of the t-norm MIN for the above given solutions + vertical accelerations 


is presented in Fig. 7. Among them the design variant of + water flooding onto the car deck in the place of broken 
L/B = 6.75 and B/d= 3.60 appears the best one. The optimum side plating 
hull form of the ro-ro ferry is shown in Fig. 8. + additional ship resistance in waves, 


has been presented. 


O The elaborated method is partly based on exact numerical 
calculations and partly on approximations obtained with 
the use of artificial neural networks. To assess the design 
merits the method based on deterministic scenarios was 
applied , and decision as to the choice of the best solution 
was based on fuzzy logic elements. It made it possible to 
present design merits in the form of liguistic variables and 
to find a design variant being a compromise between partial 
targets. 


O The optimum hull form of the ferry, resulting from 
the approximation , was verified by means of exact 
numerical methods. As results from Tab. 7 the elaborated 


5 0.350-0.360 saunas approximations appear very accurate, that can speak for that 
aN ERPE m0300.0.310 the design in question posseses appropriate (compromise) 
@0.320-0.330 =0.290-0.300 merits with regard to all partial criteria. 

Fig. 7. Values of the t-norm MIN for the best design solutions, O The presented approximations may be used to design 


CB = 0.64, CWL = 0.801, L/B = 6.72 + 6.75, B/d = 3.57+3.80 Me sce 7 
analyses (e.g. acc. [1]) or as objective functions in other 


optimization methods of ship design parameters. 


O The presented approach, due to the fuzzy logic elements 
applied to the assessing of design variant merits, makes 
it possible to simplify a multi-criterial design problem 
to that one-criterial. In such case several partial criteria 
can be represented by a single criterion , e.g. the ¢— norm 
minimizing operator, and solutions for which it reaches 
maximum can be searched for. Such approach may be used 
a.o. in optimization methods based on genetic algorithms. 
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How to divide the global resistance change 
into the components related to the pushing 
propeller and the pulling propeller 
of the double-ended ship 


Henryk Jarzyna, Prof. 


ABSTRACT 


The division of the global resistance change AR 
pushing (SR) and pulling(SD) screw propellers into the components AR, and AR,, can be 
done by taking into account two different procedures and different hypothetical assumptions 

given in the form of the division coefficients (A), and (A) „ the same in each procedure. 
In the procedure nol: AR, = AR 


In the procedure no 2: AR, = T, —(R,— Fp) (A)r and AR =T, —(R,-F,) (A), 


aap due to the simultaneous action of the 


(A) and AR, = AR,.,(A\)p, 


R+D 


Two necessary conditions were formulated to be fulfilled when selection of the coefficients (A ) Rand (A ) A 
is done. The results of the analysis are univocal. Among the possible division coefficients only one pair of 
them fulfills the formulated conditions. This pair is: 


Tr 
(ADk= m and (A 


Trap 


pea 
oe Trip 


Keywords: self propulsion tests, double-ended ships, resistance change due to propeller action. 


INTRODUCTION 


In model self propulsion tests, among other problems, 
the model resistance changes due to propeller action can 
be investigated. In the case of a ship with the pushing 
screw propeller (SR) the procedure of determination of the 
ship hull resistance change is well known from the ITTC 
Recommendation: 


AR =T-(R,-F,) (1) 


where: 

AR - ship hull resistance change 

R, - ship hull resistance when towed (from resistance test) 
T - screw propeller thrust (measured) 


F, - additional towing force (calculated) 


In special cases, e.g. double ended ferry fitted with pushing, 
SR, and pulling, SD, screw propellers being in simultaneous 
action the determination of the global resistance change can be 
performed according to the known ITTC Recommendations: 


AR pn = TR +T, (R, Fp) = Tkp (R, F,) 2) 
where: 
Tẹ Tp - thrust of SR and SD screw propellers (measured 
values) 
R, - ship hull resistance when towed (from resistance 
test) 
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F - additional towing force (calculated) 
AR,.5 - global ship resistance change due to simultaneous 
action of SR and SD propellers. 


The mentioned ITTC Recommendations do not give any 
instruction related to the necessary procedure of dividing the 
global resistance increase AR, pinto the components AR, and 
AR,, due to the both screw propellers SR and SD. Additional 
hypothetical assumptions are to be proposed to form this 
dividing procedure. 

In this paper an analysis of the possible hypothetical 
assumptions is given from the point of view of the necessary 
conditions formulated by the author. 


POSSIBLE PROCEDURES 
AND HYPOTHESES 

The division of the global resistance change AR, due 
to the simultaneous action of the pushing screw propeller SR 
and the pulling screw propeller SD into the components AR, 
and AR, can be done by taking into account two different 
procedures and different hypothetical assumptions, the same 
in each procedure. 

Two possible procedures are evident from the formal 
twofold notation of the division way. 

The first procedure (i = 1) makes direct use of the 
components AR, and AR, which are combined with the division 
coefficients (A,), and (A,),, in the form: 


AR, = ARy. (Ap 


AR, = AR (Ap 
The second procedure (i = 2) refers to the notation of global 
resistance change in the form (2): 
AR pp =T +D R, T Fp) 
and makes use of the similar notation related to the components 
AR, and AR, : 


R+D 


(3) 


R+D 


AR = Ti- (R Fok 
AR 1 Ath). 


where (R,-F,), and (R,-F,,), are connected with the same 
division coefficients (A)r and (Ab in the form: 


(4) 


(R, E Fide z (R, FCA) 65) 
(R, = Fe = (R, > FA), 
giving: 
AR, = T, —(R,—F, (A). 
(6) 


AR, =T,- R, -F (Ap 


The form of the division coefficients (A)r and (A); 
limited to a fraction of the component thrust Te slp or iie 
component P p Ppp and global thrust or power. All measured 


values (T TT T Dns Pops PoP ) are 


5 R? D, R+D’ ~ R+0? 0+D? ~ DR? DD? ~ D(R+D)? Pima. D(0+D) 
gained from the three possible types of self propulsion tests: 


¢ TestNo1 - the basic test with SR and SD screw propellers 


in simultaneous action, 
where: Tp TaT P.,P..P 


pt rip | pe! pp — measured values. 


D(R+D) 


¢ Test No2 - the auxiliary test with SR screw propeller 


only, 
where: Tgp Poe) ~ measured values. 
¢ Test No3 - the auxiliary test with SD screw propeller 
only, 
where: T p P - measured values. 


0+D? ~ D(0+D) 


The different analyzed hypotheses in the form of the division 
coefficients (A.), and (Ap are given in Tab.1, distinguished by 
the indices j, (ftom j = =1 to j= = 6). 


Tab. 1. Possible division coefficients (A), and (A), 


j = 2 thrust ratio 
To R 
j =3 power ratio 
Po Pp 


Por Pop 
j = 4 power ratio 
Polos R) Pops 0) 


=5 thrust ratio 


THE NECESSARY REQUIREMENTS 
TO BE MET 


From the division coefficients given in Tab. | these are to 
be selected which satisfy two necessary conditions formulated 
by the author. 

Condition No 1: 

The sum o the components AR, and AR, calculated by help 
of the division coefficients Ae and (Ab is to be equal to the 
measured global value AR,, 


AR, + AR, = AR. 45 (7) 
Condition No 2: 
The component AR, and AR, calculated by help of the of the 


division coefficients (A), and (A), - in each of the procedures 
(i=1 and i=2) - are to be respectively equal to: 


AR, = AR 5 
AR, = AR,, 


Theorem No 1: 

The division coefficients (A)r and (A.),, fulfill the necessary 
condition No 1 then and only ‘then when the sum of these 
coefficients is equal to one: 


(A), t (A) =] 0) 
Proof: Condition No | realized in procedure No | takes 
the from: 
AR pg, + AR p = AR Rap (A; ) pt ARR; = = 


= AR papl(A,p+ (Ajo a 


The right hand side of this notation is equal to: AR, 
then if (A), + (ADs, J=1 

Condition No Í realized in procedure No 2 takes the 
form: 


ARaot AR p= Tp- (Ro 


only 


-FP)+t Tp- (Ro- Fp) p= 


=Ta—(Ro- Fp )Aj)p+ To- Ro- FA) = 0D 
= TR+ Tp- Ro- Fp MCA jpt (A pl 
The right hand side of this notation is equalto AR, x 


then if: se +A ]= 1 because T, + Tp- R,- Fp) = Th 
~ (R, = AR uo 
Tissen No2: 


The division coefficients (A), and (A, )p fulfill the necessary 
condition No 2 then and only then if: 


T 
An and (A) = (12) 
A Trip ua Tes D 
Proof: If AR, is calculated in each of the procedures (i = 1 
and 1 = 2) then one has: 
* in procedure No 1: 


ARa =AR,, wade (13) 
* in procedure No 2: 
AR = T= (R, = Fe (14) 


and if both the values are equalized according to the condition 
No 2 then one receives: 


AR, (Ade = Tko (Ro FMA), 05 
[AR,.. + ae FIA) =T (16) 
TrolA h =T (17) 
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(18) 


If AR, is calculated in each of the procedures (i = 1, i = 2) 
then one has: 
e in procedure No 1: 


ARa F AR. wAj)y (19) 
e in procedure No 2: 
AR. = Ti (R, = FMA)p (20) 


and if both the values are equalized according to the condition 
No 2 then one receives: 


AR, o(A)p=Tp-(Ry—FyMA)p 2D) 
[AR o t R -FAD =T, 22) 
T 
(A) =- 23 
VD Tep ( ) 


The necessary condition No 2 with the demand that results 
of the division should be the same in both the procedures, leads 
to the statement that the division coefficient must be univocally 
defined as follows: 


T T 
A=- ; a=- (24) 
R Trip E Trip 


Both the coefficients fulfill simultaneously the condition 
No 1. 


CONCLUSIONS 


O Among the division coefficients given in Tab. 1, only one 


pair: T T 
Ak Ap 
i Trip j Trip 


(25) 


fulfills unconditionally both the necessary conditions: 
condition No 1: 


(A A) =1 (26) 
condition No 2: AR, ,, = AR, 
AR = AR, (27) 


O The division coefficients (A3), and (A3),,, given in Tab. 1 
for j = 3, and related to the power: 


Ppr ; (Ay) Ppp 


(A3),= (28) 
x Pore) Proves) 
fulfill the condition no 1: 
(A) t (A) = J (29) 
but they do not fulfill the condition No 2 in general: 
AR ARQ; 
(30) 
AR $ AR ps 


unless the efficiencies of the screw propellers SR and SD 
are equal to each other because in this special case only the 
following is true: 
P T P T 
DR _—- AR and DD _ _*D 
PoreD) 


(31) 


Poro) Trip Trap 
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O The simplicity of the first procedure makes it generally 
preferred in practical application. 


NOMENCLATURE 


the additional towing force 

— the parts of F, being under the influence of the SR and 
SD propeller, respectively 

— the total delivered power 

— the power of the SR propeller only being in action 

— the power of the SD propeller only being in action 

the power of the SR and SD propeller, respectively, 

being in simultaneous action 


D(0+R) 
D(D+0) 
DR? DD 


P +P 


_ _ D@+R)  D@+0) 


the algebraic mean value: P,,,= 

— the hull resistance when propelled 

— the hull resistance when towed 

the hull resistance when propelled by the SR and SD 

propellers 

— the hull resistance when propelled by the SR propeller 

only 

— the hull resistance when propelled by the SD propeller 
only 

— the parts of R, being under the influence of the SR and 

SD propeller, respectively 

the resistance increase due to the propeller action 

the resistance increase due to the simultaneous action 

of the SR and SD propeller 

the resistance increase due to the SR propeller only 

being in action 

the resistance increase due to the SD propeller only 

being in action 

— the resistance increase due to the SR and SD propeller, 

respectively, being in simultaneous action 

the pulling screw propeller 

the pushing screw propeller 

the thrust of the SD and SR propeller being in 

simultaneous action 

— the thrust of the SR propeller only being in action 

— the thrust of the SD propeller only being in action 

the thrust of the SR and SD propeller, respectively, 

being in simultaneous action 


the algebraic mean value: T „= 
Acronyms 

ITTC — International Towing Tank Conference. 
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Application of artificial neural networks to 
assessment of ship manoeuvrability qualities 


Tomasz Abramowski, Ph. D. 
Szczecin University of Technology 


ABSTRACT 


This paper presents an attempt to applying neural networks for assessment of parameters 
of standard manoeuvrability tests, i.e. circulation test and zig-zag test. Methodological 
approach to application of neural networks as well as applied network structures and neuron 
activation functions are generally presented. Also, results of simulations performed by means 
of the elaborated networks are given in comparison with test cases selected at random. In 
order to analyze and reveal general trends, correlation relationships between results from 
network simulations and test cases were calculated and are presented as well. 


Keywords: neural networks, ship manoeuvrability qualities 


INTRODUCTION 


One of the elements subjected to hydromechanical analysis 
during ship design are its manoeuvrability qualities. In 2002 
IMO adopted its MSC.137(76) resolution on “Standards for 
ship manoeuvrability” which imposed relevant requirements on 
manoeuvrability parameters of new ships. Therefore it seems 
reasonable to elaborate tools for assessing such qualities. In 
this work a method based on application of neural networks 
was proposed. Such networks have been already used for 
similar purposes. Clausen et al. [3] proposed to use neural 
networks to determine ship main dimensions in preliminary 
design stage. Koushan [7, 9] presented possible approximations 
of propeller-generated pressures and rudder forces, by using 
neural networks. Cepowski [2] applied the method in question 
to assessment of ship sea-keeping qualities. 

In this work were analyzed possible application of neural 
networks to two typical manoeuvres performed during sea 
trials of single propeller ship, namely : the circulation test 
with ship rudder laid right by 35° and the 10°/10° zig-zag test 
. The choice of the tests was conditioned by the fact of their 
standardization in IMO regulations. The presented work was 
aimed at elaboration of a method making it possible to perform 
multi-variant, simultaneous analysis of design solutions without 
necessity of multiple repetition of calculations by means of 
a simulator. 

To simulate and form the networks the MATLAB 
environment with application of the Neural Network Toolbox 


package was used [9]. A dynamic simulator being a part of 
Tribon Initial Design software was used as a source of learning 
data. The Tribon System is based on numerical solving the 
well-known differential equations of two-dimensional motion 
of ship: 


(m+m, )(i- vo -X,o7)=X, +X, +X, +X, 
(m+m, )(¥-u0-X,6)=Y,+¥e+¥, © 
(I, + J,)o +mX,(v+um)=N,+Nzp +N, 


The right-hand sides of the equations represent the inertial 
terms including added masses, and their left-hand sides — the 
hydrodynamic forces : acting on ship hull — (X Yp Np), on 
propeller — (X,), on rudder — (X,,Y,,N,), as well as due to 
such disturbances as wind and wave —(X,,, Y,,, N,.). Despite the 
learning data were supplemented by available sea trial results, 
their prevailing part was provided by the simulator.Hence 
the elaborated networks can be considered to be a simplified 
algebraic form of the Eqs. (1). 

The full set of network learning and testing data comprised 
260 cases having systematically changeable parameters. The 
set was divided in random into the learning set which took part 
in the process of calculation of network weighing factors and 
which comprised 230 cases, and the testing set which served 
solely to assess network quality and which did not take part 
in learning process and comprised 30 cases. The analysis of 
parameters in question covered the data presented in Tab. 1. 
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Tab.1. Variability range of parameters of analyzed ships 


Maximum value | Minimum value 


Parameter 


The propeller pitch ratio P/D was equal to 1 for all learning 
cases. The learning data relate to manoeuvres executed in water 
of unlimited depth. For both the manoeuvres the same rate of 
laying the rudder, equal to 2°/s, was assumed. 

The input data were normalized within the interval (1,-1) 
in accordance with the following relationship: 


(¥ nex — Yin)” & —% win) 
y = max min min (2) 
(X max — X min ) + Y min 
where: 
Yinax? Ymi © @8Sumed limits of normalization interval, in the 
considered case: y =l,y.. =-l 
n Š E max n min 
X in - minimum value in a given row of data, 
X ax - maximum value in a given row of data. 


For all elaborated networks one neuron activation function 
was applied. It has the following form: 


2 


tan sig(x; ) = E 


(3) 


where: 


x. - input value. 
1 


The run of the applied activation function is shown in Fig. 1. 


Fig. 1. The applied neuron activation (saturation) function. 


ARTIFICIAL NEURAL NETWORKS 


There is no commonly accepted, unambiguous definition 
of artificial neural network (SSN). It can be said that the SSNs 
are an attempt to connect human brain possibilities of parallel 
processing with speed and precision of computer processors. 
Such connection provides, to the networks, capability of 
learning, memorizing, calling and comparing the nonolinear, 
multi-dimensional structures of data in order to interpolate 
and extrapolate them appropriately. The traditionally applied 
modelling methods consist in collecting data and attempting 
to combine them into known physical, chemical or mechanical 
relationships. In such case to linearize mathematical models 
is necessary to make them simpler for analyzing. The method 
based on neural networks consists in an attempt to understand 
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a modeled system as an entity of relations occurring within its 
structure and it searches all linear and nonlinear relations which 
can influence its operation. Some neural networks constitute 
models of biological neural structures,but the other ones do 
not. However historically, most inspirations in the field of 
neural networks come from a desire of creating an artificial 
system capable of executing complex, possibly ‘intelligent’ 
calculations similar to those regularly done by human brain. 
Neural network structure is comprised of several processing 
units capable of mutual communicating by means of connections 
having different weighing factors. Generally, not every neural 
network has its structure similar to that shown in Fig. 2., which 
consists of the following elements: 
* Input layer — used for receiving and processing external 
signals, 
* Hidden layers (or only one layer)- whose signals remain 
within the network, 
* Output layer — used for sending signals outside. 


Each layer can be composed of several single neurons, and 
number of neurons in input and output layers is determined by 
a kind of task to be solved (number of variables and results). 

Hidden layer 


FLA 


Input layer 
moa 


Output layer 


Constant value 


C 
Constant value 


Fig. 2. Example structure of neural network 


Neurons ofparticular layers can be additionally connected with 
external elements which serve as constant values and improve 
effectiveness of network learning process. In order to clearly 
present the operational principle of neural networks, making 
use of idea of a single neuron is most effective (Fig. 3): 


xl 
Wil 
Wiz ui 
WiN 
Wio 
xn 


Constant value 


Fig. 3. Single neuron. 


Calculation of the output value y, for a single neuron is 
performed in accordance with the following relationship: 


N 
y, =F S/W, °x; + Wi (4) 
j=l 
where: 


F - applied activation function 


W, - weighing factors for particular connections 
x, - input value 
Wp - constant values attributed to particular connections. 


For given constant values and weighing factors it is easy 
to generalize the above presented relationship for an arbitrary 
number of layers and number of neurons in a given layer by 
using linear algebraic equations. The basic task in applying 
neural networks is to calculate constant values and weighing 
factors for a given network structure. The task is realized during 
learning process. The often used learning method is that called 
back-propagation one, applied in this work and presented in 
detail by Osowski [11]. 

The first operation in the learning process is to prepare 
two sequencies of data : learning one and verifying one. The 
learning sequence constitutes a set of such data which relatively 
exactly characterize a given problem. A single portion of the 
data is called the learning vector. It contains the input vector, 
i.e. the input data delivered to network input points, and output 
vector, i.e. the data expected to be generated by the network at 
its output points. After processing the input vector, the teacher 
compares the obtained values and expected ones and informs 
the network on whether the response is correct or not, if not 
— on how large error of the response resulted. The error is then 
propagated within the network but in the opposite sequence 
than that of the input vector (i.e. from output layer to input one) 
and on its basis an appropriate correction of weighing factors is 
introduced to each neuron in order to obtain a lower response 
error as a result of the repeated processing of the same input 
vetor. Such procedure is repeated again and again until the error 
generated by the network is lower than that assumed. Then the 
successive input vector is delivered to the network input points 
and the described operations are repeated. After processing the 
entire learning sequence (called the epoch) the error for the 
epoch is calculated and the entire cycle is repeated until the error 
drops under its allowable value. As mentioned above, the SSNs 
show certain tolerance for discontinuities, random disturbances 
or even small shortages in the learning set. It just results from 
their capability of generalizing the knowledge. 

The verification process of network operation follows 
the learning process. In this moment it is important to put in 
the network some specimens not coming from the training 
set, in order to check if the network is capable of effective 
generalizing the learned task. To this end is used the verifying 
sequence which has the same features as the learning one, i.e. 
its data charecterize the problem exactly and exact responses 
are known. However it is important to select the data which 
have not been used earlier for learning. This way the verifying 
sequence is presented, however the difference is that errors 
generated in the process are not propagated back but only 
number of correct responses is recorded and on this basis the 
statement is made if the network in question satisfies imposed 
requirements i.e. to which extent it has been learned. 

The initial weighing factors with which a given network 
starts learning, are numbers generated at random. After learning 
process it is always favourable to repeat all the procedure 
beginning from generation of initial weighing factors, in order 
to check the obtained results. 

In the case of large networks and learning sequences 
consisted of many thousands of learning vectors, number 
of calculations to be executed during all the learning cycle 
is gigantic, hence time-consuming. Also, it does not happen 
that a network is correctly built at once but always it becomes 
a result of many trials and errors. Moreover it can be never 
quarranteed that even a correct network does not come toa local 
minimum instead to continue finding global one. Therefore 


algorithms realizing the SSNs are fitted with mechanisms which 
make it possible to control speed and quality of learning. These 
are the so called learning and momentum coefficients. They 
influence steepness of activation function and control speed of 
influence of change of weighing factors on learning process. 


CALCULATIONS 
FOR CIRCULATION TESTS 


To simulate the circulation manoeuvre a network of 9x12x7 
structure was elaborated. The input values to the network 
(arguments) were in turn the following : the ship length 
between perpendiculars, L,,, ship breadth B, ship draught T, 
block coefficient of ship hull underwater part, C,,, longitudinal 
position of ship centre of gravity, L.,, propeller diameter D, 
ship speed at the beginning of the circulation manoeuvre, V, pp 
rudder aspect ratio À and the total rudder area A. 

The neural network outputs were the following : transverse 
and head translations (transfer and advance) at ship course 
change by 90°, tactical diameter of circulation and head 
translation for the course change by 180°, steady circulation 
diameter, total decrease of ship speed at course change by 360°, 


V/V >» a8 well as the drift angle B during steady circulation. 


The elaborated structure is presented in Fig. 4. 


Fig. 4. Neural network for circulation manoeuvre 


In Fig. 5 the run of the network learning process is shown. 
It can be observed that after about 300 epochs of learning 
the network did not exhibit any further improvement of 
approximation quality. 
10° Performance is 0.00186414 


Performance 


0 50 100 150 


20 20 30 0 


350 Epochs 


Fig. 5. Run of network learning process for circulation manoeuvre 


The simulation of the manoeuvre was performed on the basis 
ofa case selected from the testing set. The results of the simulation 
compared with the expected values are presented in Fig. 6. 

It can be observed that for the course angle change by 90° 
the network error is of the order of the ship breadth. In the 
case ofthe course angle change by 180° the error is somewhat 
smaller - equal to about a half of the ship breadth. The drift 
angle at steady circulation is of a similar order. The difference 
of diameters at steady circulation is almost imperceptible. For 
purposes of this presentation one case characterized by the 
greatest errors was selected out of 30 testing ones. 
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Head distance/Ship length 


1.0L 2.0L 
Transverse distance/Ship length 


3.0L 


Fig. 6. Results of simulation of neural network for circulation manoeuvre 


In Fig. 7. is presented the correlation between the drift 
angle values determined by means of neural networks and its 


R2= 0.9354 


N 
=] 


D 29 


Neural Network prediction 


N 


10 1l 12 13 14 15 16 17 18 19 2C 
Test case value 
Fig. 7. Correlation between values of the drift angle at steady circulation, 
obtained from the network, and relevant testing values. 


testing values, and in Fig. 8 — the correlation of values of the 
head translation at the course change by 90°. 

Example results of network prediction values compared 
with relevant testing ones are shown in Tab. 2. The values 
concerning the translations (distances) are related to the ship 
length. Apart from the tests of correlation between results 
obtained from the network and testing cases, also capability 


Tab. 2. Example simulation quality of neural network used for assessment of circulation parameters 


Parameter 


Advance/L [-] 


Course change 


Transfer/L [-] 


1.67 | 1.98 


Tactical diameter [-] 


3.58 | 4.19 


Advance/L [-] 


Testing values 


Diameter [-] 


2.32 | 2.86 


Steady circulation V/Vapp [-] 


| 
| 
2.13 | 2.34 
| 
| 


0.38 | 0.44 


Drift angle [°] 


Advance/L [-] 


15.29 | 12.96 


3.23 | 3.71 


Transfer/L [-] 


1.65 | 1.94 


Tactical diameter [-] 


3.57 | 4.12 


Advance/L [-] 


Diameter [-] 


2.34 | 2.86 


E 
=l 
£ 
D 
(= 
T 
x 
= 
E 
Z 


Steady circulation V/Vapp [-] 


| 
| 
2.09 | 2.32 
| 
0.38 | 0.46 


Drift angle [°] 


Relative error in [%] of a given value 
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15.63 | 12.95 


na 
h 


presented correlation between the test values for the 1st and 
0.9741] ¢ 2nd overshoot angle and the relevant values approximated by 
the neural network. 
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Fig. 10. Neural network for simulation of zig-zag test 
10° Performance is 0.00462386 
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Fig. 8. Correlation between values of the head translation at the course 
change by 90°, obtained from the network, and relevant testing values. 


Performance 


of the elaborated networks, related to modeling the general 
trends of circulation parameters, was analyzed. In Fig. 9 is 
presented an example of the testing of influence of the ship 
hull block coefficient C, on circulation parameters by using 
the elaborated networks. 


n w © w w 500 
504 Epoche 


Fig. 11. Run of learning process of the network for zig-zag test 


Aa 1,Aa 2 — differences 
in overshoot angles between 
neural network prediction 
and the test case 


Test case value 


10 
Neural network prediction 


STBD rudder/heading 


Port rudder/heading 
S 


Test case value 


Fig. 12. Results of zig- zag test simulation by using the neural network 


1.0L 2.0L 3.0L 4.0L 5.0L 12 


Fig. 9. Testing of the trend representing capability of the elaborated R*= 0.8910 
networks. Influence of ship hull block coefficient. 


CALCULATIONS FOR ZIG-ZAG TESTS 


_ 


© 


wo 


To simulate the zig-zag test the network of 9 x 11 x 5 
structure was elaborated. Inputs to the network were identical as 
for the circulation test. It was the following parameters: duration 
time of the course change by 10° (initial turning time), t,; Ist 
and 2nd overshoot angle, Ay, Ay, ; 1st and 2nd time to check 
yaw (time of rudder inertia angle), t,, t,. The network structure 


a 19 “9: 6 


Neural Network prediction 
oe 


is shown in Fig. 10. The run of the network learning process 
is shown in Fig. 11. The learning process was terminated after 5 
600 epochs. 
In Fig. 12 are presented results of the zig-zag test simulation 4, 5 6 = g 9 10 1 
compared with those of the test case. Errors in representing Test case value 
particular time values and overshoot angles are contained within Fig. 13. Approximation quality of the neural network, based on the 
a few percent interval of relevant values. In Fig. 13 and 14 is correlation for the Ist overshoot angle 
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Neural Network prediction 


20 


The example results of the zig-zag test prediction, compared with relevant test values are given in Tab. 3. The obtained 
prediction quality is somewhat worse than that in the case of simulation of circulation manoeuvre. 


Tab. 3. Example simulation quality of the neural network used for assessment of zig-zag test parameters 


Parameter Values 


1* overshoot angle. Ay, [°] 


2"* overshoot angle. Ay, [°] 


Initial Turning Time. t, [s] 


Test values 


1* Time to check yaw. t, [s] 


2"! Time to check yaw. t, [s] 


1* overshoot angle. Ay, [°] 


2"* overshoot angle. Ay, [°] 


Initial Turning Time. t, [s] 


1“ Time to check yaw. t, [s] 


Neural network 


2™ Time to check yaw. t, [s] 


Relative error in [%] 
of a given value: 


4 6 8 0 2 4 K 
Test case value 

Fig. 14. Approximation quality of the neural network, based on the 

correlation for the 2nd overshoot angle 


CONCLUSIONS 


Essence of the presented method consists in its possible 
applications to optimization and multi-criterial automatic 
design of ships. A difficulty in applying traditional 
manoeuvrability simulators is their low capability of 
fast and simultaneous analyzing many design solutions 
within a similar time interval. This is the feature which 
must characterize optimization algorithms. Sometimes the 
feature is more desirable than a high accuracy of results 
achieved at the expense of large time consumption. In the 
case of optimization algorithms it is necessary first of all 
to distinguish a better solution from worse one. 
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O The presented method makes it possible to simultaneously 


determine parameters of analyzed manoeuvres depending 
on design parameters of ship. The presented correlation 
diagrams show that the networks correctly recognize trends 
and that the obtained results are not accidental. 


O The simulation quality for circulation manoeuvre is a little 


better than that for zig-zag test, though it would be possible 
to improve the simulation quality by supplementing the 
learning set. 


O It should be highlighted that the presented method is not 


aimed at improving ship motion modeling quality but 
only at improving and adjusting the existing algorithms to 
demands of design process. 


O The presented neural networks are capable of exact 


interpolating however their capability in extrapolating 
have not been sufficiently recognized so far, hence during 
application of the method any results obtained for values 
of the ship parameters being beyond their intervals given 
in Tab. 1, should be considered very cautiously. 


O However from another side the traditional manoeuvrability 


simulators make as a rule use of regression relationships to 
determine force coefficients. The method in question makes 
it possible to continue research in the presented direction, 
e.g. by taking into account corrections associated with water 
depth, ship trim or application of multi-propeller propulsion 
systems. 


O It is also possible to use neural network for dynamic 


simulation of manoeuvre; in this case for a given ship as 


network inputs would be used its kinematic parameters 
at a given time instant t together with taking into acount 
instantaneous position of rudder and screw propeller setting, 
and network outputs — kinematic parameters after a given 
time interval At. 


O Certainly further investigations on the method would be 
carried out in the direction of bettering the learning data 
by supplementing them by results from real sea trials, and 
with a view of presenting example optimization results. 
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for multiaxial stress in viscoelastic solids 


Janusz Kolenda, Prof. 
Naval Academy of Gdynia 


ABSTRACT 


On the basis of Hooke law for multiaxial stress in elastic solids, similar relationships for 
viscoelastic materials are considered. It is assumed that the material is homogeneous and 
isotropic, and that the Kelvin-Voigt’s model is applicable to normal strain components. An 
analogous model is also taken for shear strain components. It is shown that the ratio of 
coefficients of viscous damping of normal and shear strain components is equal to the ratio 
of Young modulus and shear modulus. As a result, the modified Hooke 5 law for multiaxial 
stress in viscoelastic materials has been formulated which includes three material constants: 


Young modulus, Poisson 5 ratio and coefficient of viscous damping of normal strain. 


Keywords: Hooke’s law, viscoelastic material, Kelvin-Voigt’s model, multiaxial loading 


INTRODUCTION 


For the purpose of solving physical problems considered 
in the present paper, it is assumed that an engineering detail is 
made of a homogeneous isotropic metal and loaded below the 
yield point. The elastic behaviour of a particular homogeneous, 
isotropic metal at a given temperature is completely defined 
by the Young modulus and the Poisson’s ratio [1, 2]. However, 
even in the region below the proportional limit metals are 
not perfectly elastic. If a load is suddenly applied and then 
maintained constant, a small amount of “creep” will always 
be observed [1]. This creep is due to anelastic strain which 
may arise from any of several sources, such as the presence 
of grain boundaries, twin boundaries, or slip bands, the 
diffusion of interstitial solute atoms and the phenomenon of 
magnetostriction. Particular attention is focused on anelastic 
strain in the vibration theory of continuous systems [3, 4]. To 
exemplify the problem, let us mention that contrary to free 
vibrations of perfectly elastic systems where the principle 
of energy conservation may be assumed, in reality the initial 
energy of isolated vibrating systems diminishes with time as 
a result of energy dissipation caused by damping properties 
of structural materials. This phenomenon is very complex but 
there are several approximate models of anelastic materials and 
damping mechanisms which are simple enough to be accepted 
in practice, for instance the Kelvin-Voigt's model (spring and 
dashpot in parallel) and hysteretic damping [3-6]. In the present 
paper, the Kelvin-Voigt’s model is applied to viscoelastic 
materials under three-dimensional loads. The second part of 
this paper, to be published separately, will be devoted to the 
applications of relationships contained herein, with an emphasis 
made on the behaviour of viscoelastic solids under multiaxial 
loads and dissipation energy. 
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HOOKE’S LAW 


When structures, machines, and engineering details are in 
service, they are usually subjected to surface forces or body 
forces (inertial, gravitational, or electromagnetic) that cause 
combined stresses in their elements. To completely define an 
element, it is necessary to specify the components of the stress 
tensor represented by the array: 


Ox Tx Tax 
T=|t, O, Ty (1) 
| Txz Toz O, | 


where o are the normal stress components and z are the shear 
stress components on three orthogonal planes passing through 
the point. If surface couples and body couples may be ignored, 
the number of independent stress components reduces from nine 
to six since then t =T,,1, =t , andt =t_ Consequently, 
. xy yx Z zy’ zx XZ 
deformation produced by stress components can be completely 
specified by six strain components: normal ones £, E €, and 
shear ones y., Y, „Y, The Hooke’s law states that in elastic solids 
xy? !yz? TA i è 
between these components the following relations hold [7, 8]: 


Ee.=o,= vlo, + o,) 
Es, =0, -vlo +0,) 
Ee, =0,- vlo, + o,) 

CY xy = Ty 

Oy, = TW, 

Gy, = Tax 


(2) 


Here E is the Young modulus, v is the Poisson’s ratio, 
and: 
E 


2(1+v) (3) 


is the shear modulus. 

For any combination of loads, there exist at any point three 
orthogonal planes upon which only normal stress components 
act. These are the principal stresses, designated in declining 
order of magnitude: o, 6,, and o,. Therefore, for the particular 
case where x, y and z are the principal directions, Eqs (2) 
yiel: 

Eg, =0, -v(o, + o;) 


Ee, =o, - vlo; +0) 


Es, = 6, — vlo, + o) e 


Yio = Yz = Ya = 9 


From Mohr’s stress and strain circles [7, 8] it follows that 
the maximum shear stress and the maximum shear strain are: 


t=0;5(0 —6,) (5) 
y= 26, (6) 

The relationship between t and y reads: 
Gy=t (7) 

For the given stress state that: 
6 = =6,7;5;—-0 (8) 
one gets: 

Fe, =¢,—0(0'—«,)=0,0 +0) (9) 


TO, 


(10) 
HOOKE’S LAW FOR VISCOELASTIC 
MATERIALS 


Let us consider what happens to a long rod of homogeneous 
isotropic material ifa tensile force is applied to the end of the rod 
in the direction x parallel to its length. Of course, this rod will 
exhibit a uniform tensile strain e, produced by the stress o. 


For elastic state: 
Es =o 
X X 


(11) 


Also, it will be observed that the fractional change in 


diameter is proportional to the axial strain and such that: 
E€ =— VE ,E€ =—VE 
y x Zz x 


(12) 


If the force and, consequently, the stress o, is time- 
dependent: 


o = 6 (t) (13) 
it is assumed that Eqs (12) are also valid, i.e.: 
e (t) =- ve, (t) , e(t) =— ve (t) (14) 


Differentiation of Eqs (14) with respect to time gives: 


é(t)=—vé,(t) , é,(t)=-ve,(t) 09 
In order to account for viscoelastic properties of the rod in 


this case by means of the Kelvin-Voigt’s model, Eqs (11) and 
(12) must be replaced by: 


Es, +né, =9, (16) 


=-vE, (17) 


where y is the coefficient of viscous damping of normal strain, 
and [9]: 
(18) 


In Eq. (18), v, is the Poisson’s ratio in the case of anelastic 
strain, v, is related to elastic deformations, v is caused by 
microslips, and v, reflects the influence of microcracks. If the 
load is such that the solid remains in elastic state, then: 


VA 


E€ =—VE ,E 
y a x 


v =v ty +o 
a e S c 


v, =V, =v, V 5V, 50 (19) 

Since experiments show that Eq. (19) corresponds very 
closely to the behaviour of many solids under small loads [9], 
in what follows Eqs (12), (14) and (15) will be applied. 

Now we shall compute the normal strain in the x-direction 
caused by normal stresses in the y- and z-directions, where x, 
y and z are any three mutually perpendicular directions. 

For the strain due to o, one gets analogously to Eq. (16): 


Ee, +né, =o, (20) 
Thus: 
VEE, + ons, = vO, (21) 
But in this case: 
€, =—-ve, , $, =—v6, (22) 
so that: 
Fe, + 1, =—v0, (23) 


Similarly, for the stress in the z-direction one obtains: 
Ee, +1é, =—-v0, (24) 
Eqs (16), (23) and (24) are linear. These equations and 
principle of superposition lead to the following equation for 
the strain in x-direction: 
Ee, +n, =9, — vlo, 45) (25) 
Equations for the strains in the other directions can be found 
in the same way. 
When x, y and z are the principal directions, Eq. (25) 
becomes: 


Eg, +n =9, — vlo, + o) (26) 
For the stress state described by equations: 
ot) ==0 t); 00 (27) 
we have from Eq. (26): 
Es, +n, =9, (I + v) (28) 


Now let us assume that for viscoelastic solids loaded below 
the yield point, the maximum shear stress and strain can be 
expressed by Eqs (5) and (6), determined for elastic solids. 

Differentiation of Eq. (6) with respect to time yields 

y= 28, (29) 

Analogously to Eq. (16) for normal strain component, for 
viscoelastic materials and shear strain component instead of 
Eq. (7) one can write [3, 4]: 


Gy+AV=T (30) 


where A is the coefficient of viscous damping of the shear 
strain. 
By Eq. (28): 


2Ee, + 2né, = 20,(1+v) (31) 
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Through Eqs (6), (10) and (29), Eq. (31) takes on the 


form: 


Ey +ny = 21(1+v) (32) 
that is: 
n ‘ 

+——y=T 

aren (33) 
From comparison of Eqs (30) and (33) it is seen that: 

a l r 
2(1+0) 64) 


Summarizing, the modified Hooke’s law for multiaxial 
stress in viscoelastic solids can be expressed as: 


Es, +16, =o, — vlo, +o,) 


Es, +né, =0, - vlo, +0,) 


x 


Es, +né, =9, - vlo, +o, ) 


; (35) 
Gy + AY wy = Ty 
GY,- + ce =Ty, 
Gy, tAV~. = Tx 
where: 
—_ E Nn 
al+v) > — 2(1+v) (8) 


Another derivation of Eqs (34) and (35) can be found in 
[10]. 
Eqs (35) can be rewritten in the following alternative 


forms: 
Es, +né, =6, - vlo, +6,) 
Ey, + Wy =20+v)t, (37) 
Ey,, + NVy, = 2(1+ vt, 
Ey x +1 x = 2(1 + via 
or 
= ( — v)\(Es, i né, )+ vlEle, F e,)+ nlé, oe é,)| 
E (1 + v)(I — 2v) 
a (1 T v)(I - 2v) 
(1+v)(1-2v) 
xy 2(1 + v) 
T,= Ey, +N 
yz 2(1 + V) 
= 2(1+v) 


When x, y and z are the principal directions with respects 
to stress, Eqs (37) become: 
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Ee +e =5,= vlo, + 5,) 


Es, + es =6, — vlo, +0) 


Es, + né; = 0, - vlo, + A) 
Ey, +N =0 
Ey», +NY23 =0 
Ey; +13, =0 


It means that, contrary to the stress-strain relations (4) in 
perfectly elastic solids, viscoelastic materials subjected to 
normal stress components may exhibit not only normal but 
also shear strain components (for instance after removal of 
a shear load). 

Eqs (35) and (37) through (39) present the modified 
stress-strain relations in viscoelastic solids under triaxial load 
conditions. In what follows those in simpler load cases are 
listed. 


(39) 


Plane stress (o, = i =) = 0) 


Ee, +1é, =O, — v0, 
Ee, +né, =, — v0, 
Fe, + nė, = -vlo, +0,) 
EY) +My = 21+ v0) 
BY. + NY yz = 9 
By rit, =o 


and 


a Es, + te; +v(Ee, + nė, ) 


x 


(40) 


10" 
a Bey +N, +v(Ee, +n€,) 
ee 


o -Pe tiy 
Y a(l+v) 


Plane stress (£, = y =0) 


yx 


Es, +n, = (+v), (i-v)-vo, | 
)-vo,| 
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Ey. +NY%xy =0 (44) 
Eya 02 = 9 
By, + ny,, =0 
and 
o, = Ee, +n, (45) 
Shear stress ( Tt, £0) 
Es, +né, =0 
Ee, + né, =0 
Es, +neé, =0 (46) 
Evy + May = 201+ v) ty 
Ey, + Y= 0 
Ey, + ny, =9 
and 
o -Paty (47) 
T (+o) 


As it is seen, the two-parameter Kelvin-Voigt’s model 
of viscoelastic materials leads to relatively simple stress- 
strain relations. The stress-strain relations based on the 
three-parameter model of viscoelastic materials (spring and 
Maxwell’s model in parallel) are given in [11]. 


CONCLUSIONS 


O The Hooke’s law for multiaxial stress in elastic solids has 
been modified to account for damping properties of real 
materials. 


O The Kelvin-Voigt’s model for viscoelastic materials has 
been applied to normal and shear strain components. 


O The viscoelastic behaviour of a homogeneous, isotropic 
material at a given temperature is completely defined by 
Young modulus E, Poisson’s ratio v and coefficient n of 
viscous damping of normal strain. 


O The aforementioned constants can be determined in 
a uniaxial test. 


O Itis shown that 
E 
N- Ž=2(1+v) (48) 
A G 
where G is the shear modulus and à is the coefficient of 
viscous damping of shear strain. 


NOMENCLATURE 


— Young modulus 

shear modulus 

— time 

three mutually perpendicular directions 

— shear strain component 

— normal strain component 

— coefficient of viscous damping of normal strain 
coefficient of viscous damping of shear strain 
— Poisson’s ratio 

— normal stress 

— shear stress 

principal directions 
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ABSTRACT 


Application of steel sandwich panels to ship structures requires many problems to be 
solved. Joints between the panels as well as those between the panels and other structures 
is one of the more difficult problems associated with the structures in question. This paper 
presents the searching for process of optimum geometry of a panel-to-panel joint of 
longitudinal arrangement, performed by means of the ANSYS software. A configuration 
was searched for of parameters which can ensure as-low-as possible values of geometrical 
stress concentration coefficients at acceptable mass and deformations of the structure. 


Analysis of the obtained results made it possible to propose the optimum geometry of the 
considered joint. 


Keywords: steel sandwich panels, connection, optimization 


INTRODUCTION 


Ship structures are evolving with time. The need for 
increasing ship reliability, lifetime and safety results mainly 
from economic motivations. Along with increasing prices of 
steel and fuels, which exert high impact on decisions made in 
shipbuilding industry, one of the main impulses accelerating 
development of new, better structures is a tendency to make 
the shipbuilding time shorter and, at the same time, the time 
between its inspections and repairs longer. The pressure towards 
shortening the shipbuilding time is so high that shipyards have 
been transformed in sort of assembly factory of ship sections 
produced in smaller specialised cooperating plants. The next 
factor which is becoming more and more important in decision- 
making in the shipbuilding industry is ecological safety. The 
need for protection against ecological disasters provoked by 
ship accidents (collisions, groundings) was the main motivation 
for designing two-skin hulls. The development of the already 
existing structures is mainly done by evolution - structure 
improvement based on many-years’ experience gained by the 
Classification Societies. In successive regulations higher and 
higher requirements are formulated (especially those concerning 
the ecology), with an intention to reduce the risk relating to 
ship operation. Along with the improvement of classical ship 
structure, the investigations are in progress over completely 
new structures which can be used as ship components. A list 
of these new structures includes steel sandwich panels, which 
reveal extremely favourable bending resistance and stiffness 
in relation to mass. 
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STEEL SANDWICH PANELS 


Briefly speaking, the sandwich panels are two-skin 
structures revealing favourable stiffness and bending resistance 
in relation to mass. A concept of multi-layer structures has been 
well known for years and is in common use in building of plastic 
watercraft (laminate-foam-laminate, for instance). Employing 
advantageous multi-layer structures in shipbuilding was 
stopped for decades by technological problems (deformations 
generated during welding of extremely thin plates to their 
stiffeners, among others). And only recently the development 
of laser welding technology has made it possible to produce 
steel sandwich panels (Fig. 1). 


Fig. 1. Industrial production of steel sandwich panels [4] 


The laser welding technology was introduced to the industry 
about ten years ago. The need, having its origin in economic 
motivations, to reduce the mass of the structure is so high that 
it compensates high costs of laser welding. 

But before making full use of the advantages offered by 
the steel sandwich panels, some problems concerning their 
application are to be solved. In particular, problems which stop 
their use as carrying elements of a structure concern [1, 2, 3]: 
strength 
fatigue 
fire and acoustic insulation 
transmission of vibrations 
reaction to dynamic and point loads 
corrosion resistance 
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The use of sandwich panels in ship structures also requires 
analysing issues referring to the presence of cables and 
pipelines, communication holes, repairs of damages, operating 
control and utilisation. A separate group of problems, of high 
significance for the final shape of the structure, are joints 
between the panels, and between the panels and the classical 
structure. 


SELECTED PANEL JOINTS 


Joints between panels or between a panel and the classical 
structure are difficult problems affecting practical application 
of the panels. Solving a number of minor problems concerning 
these joints is a basic condition for their future use as 
components of larger structures. There are a variety of types 
of joints, including, for instance, panel-panel or panel-classic 
structure joints, cross-joints, panel height reduction, and panel- 
support joint. One of more common constructional centres in 
the panel structure is the front panel-panel joint in longitudinal 
arrangement, in which the welds which link the panels are 
parallel to the direction of the stiffenings, see Fig. 2. 


Panel joint area 


Fig. 2. Positions of I-core panels in longitudinal arrangement. 


Numerous ways of execution of these joints are proposed. 
Sample (but not only possible) solutions are shown in Fig. 3. 
Individual proposals of these joints differ by mass, expected 
stiffness, cost, and time of execution. Within the framework 
of preliminary studies, two panel joints were selected for 
further examination from the point of view of technological 
realisability. In the assessment a number of factors were taken 
into account including: number of welds, preparing edges 
for welding, width of groove, assembly tolerance, ability to 
cut edges after fitting, fixing elements during the assembly, 
access to places be welded, smoothness of the joint, ability 
to prefabricate, and expected after-welding deformations. As 
a result, the joints which turned out most favourable from 
the point of view of technological realisability were those 
making use of cover plates and rectangular profiles, see Figs 
3c and 3b. 


a) b) 
c) d) 


e) 2] 


Fig. 3. Proposed joints 


SEARCH FOR OPTIMUM GEOMETRY 
OF THE SELECTED JOINTS 


To assess strength-to-mass ratios of the selected joint 
geometries, reduced stresses were calculated, using the Finite 
Element Method (FEM), in the joints loaded with stretching 
and compressing forces acting in the direction perpendicular 
to the joint. 

The geometry of the selected joints was described using 
parameters shown in Fig. 4, in which the adopted coordinate 
system and the symmetry plane are marked as well. It was 
assumed that the joined panels have the same geometry. 
Because of repeatability of the joint structures, the “distance 
from panel end to first stiffener” parameter takes the value of 
4 or ⁄2 of the distance between the stiffeners. The parametric 
description was formulated in order to search for such 
a configuration of parameters for which the joint reveals the 
lowest coefficient of geometric concentration at relatively 
small mass and deformation. The reduced Mises stresses 
were used as the measure of joint’s ability to carry the load 
(stretching or compressing force). The preferential criterion 
of joint evaluation was the stress, at an acceptable mass of the 
joint (up to 10% of panel mass). If, for some configurations 
of parameters, joints revealing similar stresses were obtained, 
the joint with smaller mass was selected. Two load conditions 
were assumed, which were stretching and compression in the 
panel plane. 


NUMERICAL MODELS OF THE JOINTS 


The calculations were performed in ANSYS environment. 
The assumed symmetry of the geometry and load, Fig. 5, has 
made it possible to analyse 1⁄4 of the model, which considerably 
speeded up the calculations. To compensate possible effect of 
panel deflection on the operation of the joint, loads were applied 
at a distance equal to 2.5 distances between the stiffeners. 
Writing the geometry in the parametric form provided 
opportunities for performing a large series of calculations with 
changing combinations of parameters, without the necessity to 
build a new model for each combination. The modelling was 
done using the PLANE183 element- a plane solid-type element 
with the square shape function. The models had a carefully 
selected regular grid consisting of about 15 thousand elements. 
The side length of an individual element was approximately 
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Fig. 4. Parameters describing the geometry of selected joints 


equal to 0.1mm. The steel material constants were assumed for 
the isotropic model: E = 2e° MPa and v = 0.3. 


symmetry axes 


Fig. 5. Panel loaded with stretching force (P) and compressing force (-P). 
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Fig. 6. Calculation model with assumed boundary conditions. 


A sample model of the joint is given in Fig. 6, along with 
the applied boundary conditions and load. During numerical 
modelling of sandwich type structures, the most difficult and 
the most questionable model elements are laser welds. The 
technology of welding the stiffeners to the plating assumes laser 
welding from outside. Specific nature of the laser-welded joints 
is the reason why some gaps form between the stiffeners and the 
plating (their dimension in the discussed panel structures can 
approximately reach 0.1mm). In the ending areas of the welds 
technological notches are localised, which, when geometrically 
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modelled in the numerical model, would lead to unrealistic, 
excessively high concentration of stresses [5, 6, 7, 8, 9]. For the 
computer model to reflect the real laser joints in a most realistic 
way, so-called circular concentrators are introduced to replace 
the welds in those places, see Fig. 7. These concentrators create 
geometric smoothening of the notches in the form of circles 
with specially selected dimensions and positions [8, 10, 11, 12]. 
Of extremely high importance in weld modelling is taking into 
account different nature of laser weld operation compared to 
a joint made using a traditional technology. The laser weld are 
characteristic of three operation stages [1]: 1 — rotation of the 
joined elements; 2 — contact of the plating with the stiffener; 3 
—common displacement of the plating and the stiffener in the 
same direction, but with the load carried by a much smaller 
cross-section than in classically welded structures. In order to 
control the nature of weld operation, relevantly defined “contact 
regions” have been introduced at selected stiffener edges. 


Fig. 7. Laser welded joint [1] and the FEM model. 


RESULTS 


To make the interpretation easier, the results of calculations 
are presented in the form of diagrams. In these diagrams the 
vertical axis represents the geometric concentration coefficient 
Kg, calculated as the ratio of the maximum of the reduced Mises 
stress to the nominal stress. On the horizontal axes, different 
variable parameters describing the geometry of the joins are 
shown. Each configuration of parameters is represented by 
one point in the diagram. A huge volume of results has made 
it possible to build a surface which illustrates the effect of 
particular parameters on stress concentration. These surfaces 
were smoothed using the distance-weighted smallest square 
method, for which the calculations were done by the code 
“Statistica”. The data discussed in the article solely refer to 
stresses. 

Fig. 8 presents the results concerning the joint making use 
of cover plates. When increasing the “half cover plate width” 
parameter the stress concentration decreases and reaches a local 
minimum for 20-30mm. Further increasing of the plate width 
results in stress increase by about 40% in the vicinity of the 
first stiffener. Within the parameter range of 10-40 mm the 


deformations in the y-axis direction increase for stretching 
load and decrease for compression. When the width of the 
cover plate exceeds the distance to the first stiffener (40- 
50 mm), the Mises stresses rapidly decrease. Deformations in 
the y-axis direction within this range decrease for the cover 
plate thickness of 2.5 and 3 mm, and reach a local minimum 
for 50 mm, but increase for the cover plate thickness of 2 mm. 
More detailed analysis of the effect of particular parameters 
on stress and deformation of the joint has made it possible to 
select the most favourable configuration of parameters. The 
optimal configuration of parameters for the panel joints having 
the “distance from panel end to first stiffener” parameter equal 
to 4 turned out to be the joint with the cover plate thickness 


equal to 2.5 mm and the parameter “half cover plate width” 
equal to 50 mm. 


Diagram of geometric concentration stress coefficient Kg 
for cover plate panel joint with the distance form panel end 
to first stiffener equal to 1⁄4 


cover plate sail a$ ae 
thickness [mm] ~ „S half cover plate 
D> oy width [mm] 
Diagram of geometric concentration stress coefficient Kg 


for cover par panel joint with the distance form panel end 
to first stiffener equal to 14 


cover plate 2” A “i 
thickness [mm] 7 œ __ halfcover plate 
a® «0 width [mm] 


Fig. 8. Geometric concentration coefficient as a function of geometric 
parameters of the cover plate joint - stretching 


For the cover plate panel joints with the “distance from 


panel end to first stiffener” parameter equal to 1⁄2, the maximal 
values of the reduced Mises stresses their reach minima in the 
vicinity of the parameter value equal to 40 mm for all cover 
plate thickness cases when changing the “half cover plate 
width” parameter. Largest deformations are observed for the 
parameter equal to 30 mm when compressing. When compared 
to the 4 case, panel joints with the “distance from panel end to 


first stiffener” parameter equal to 4 reveal five times as large 
deformations in the y-axis direction during compression, at 
a comparable level of the maximum Mises stresses. The optimal 
configuration of parameters of the analysed joint is the cover 
plate joint with cover plate thickness equal to 3 mm and the 
“half coven plate width” parameter equal to 40 mm. 


Diagram of geometric concentration stress coefficient Kg 
for rectangular profile panel joint with the distance 
form panel end to first stiffener equal to 1⁄4 


Diagram of geometric concentration stress coefficient Kg for cover 
plate panel joint with the distance form panel end to first stiffener 
equal to 1⁄4. Kg[-] = distance weighted smallest square smoothing 


t © 
æ% mass of joint [%] 


Fig. 9. Geometric concentration coefficient as a function of geometric 
parameters and mass of the rectangular profile joint - stretching 


For rectangular profile panel joints with the “distance from 
panel end to first stiffener” parameter equal to '4, the maximal 
values of the reduced Mises stresses and deformations in the 
y-axis direction depend much less on whether the joint is 
stretched or compresses than it was observed for the cover 
plate joints, Fig. 9. For small values of parameter p3 (profile 
insertion) the difference between deformations in the x-axis 
and y-axis directions are visibly smaller that for the cover plate 
joints, but when the profile insertion increases these differences 
also increase (still remaining much smaller than for the cover 
plate joints). After comparing stresses and corresponding 
masses of particular joints we can see that for the same mass the 
maximum reduced Mises stresses change nearly by twice when 
changing the parameters. Also at the same stress level, changing 
parameters can result in mass change by as much as 20%. The 
optimum configuration of parameters is a joint with the gap 
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x = 3.25 mm, profile thickness equal to 3 mm and profile 
insertion equal to 5 mm. This joint reveals the lower stress 
level at the mass equal to 9% of panel mass. Moreover, the 


at comparable masses of the both joints (within the ranges of 
their optimal parameters) the rectangular profile joint secures 
much lower stresses and deformations, in particular in the 


deformations for these parameters also are the smallest. 


Diagram of geometric concentration stress coefficient Kg 
for rectangular profile panel joint with the distance form panel end 
to first stiffener equal to 34. Distance weighted smallest square smoothing 


Diagram of geometric concentration stress coefficient Kg for cover plate 
panel joint with the distance form panel end to first stiffener equal to 1⁄4. 
Distance weighted smallest square smoothing 


Fig. 10. Geometric concentration coefficient as a function of geometric 
parameters and mass of the rectangular profile joint - stretching 


The results obtained for the rectangular profile panel joints 
with the “distance from panel end to first stiffener” parameter 
equal to 2 are given in Fig. 10. The optimal configuration of 
parameters is the joint with the gap x = 5 mm between the joined 
panels, profile thickness equal to 2.5 mm and profile insertion 
equal to 5 mm. This joint reveals the lower stress level at the 
mass of joint equal to 8% of the panel mass. 

The analysis of the above results makes it possible to 
formulate some summarising conclusions and compare the two 
examined types of joint: the cover plate joint and the rectangular 
profile joint, Fig. 11. 

The cover plate panel joint reveals much higher geometric 
concentration coefficients. Although it has solution intervals 
characterised by much smaller mass, maximum stresses in those 
intervals reach several times as high levels as those observed in 
the rectangular profile joints. Also the maximum deformations 
are much smaller for the rectangular profile joints. To sum up, 
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y-axis direction. 


Comparison diagram of cover plate and rectangular profile joints. 
Panel joints with the distance from panel end to first stiffener equal 
to 1/4 for stretching. Distance weighted smallest square smoothing 


p3 [mm] 
10 x cover 


half cover plate width 


Comparison diagram of cover plate and rectangular profile joints. 
Panel joints with the distance from panel end to first stiffener equal 
to 1/2 for stretching. Distance weighted smallest square smoothing 


p3 [mm] 
10 x cover am 
plate thickness A x [mm] 
half cover plate width 


Fig. 11. Comparing the cover plate joint (points) 
and the rectangular profile joint (surface) 


The final solution of the panel-panel heading joints in 


longitudinal arrangement is the rectangular profile joint shown 
in Fig. 12. 


Oe 


Fig. 12. Patterns of reduced Mises stresses and deformations 
(in 600:1 scale) for the optimum design (compression). 


SUMMARY AND CONCLUSIONS 


O The development of ship structures, dictated by strong 


economic factors, leads to better and better solutions. 
At the same time, concepts are formulated to replace 
fragments of classical ship structure by new-type elements. 
New solutions which can lead to remarkable reduction in 
mass and production time of the ship are steel sandwich 
panels. At present they cannot be used as ship bearers, 
because of the lack of relevant regulations of Classification 
Societies. But the economic pressure is so high that it 
forces increasing the use of panels in ship constructions. 
The panels are used for building between-decks, systems 
of ramps, decks of river vessels, balconies on passenger 
ships, stairways, and/or partition walls. Shipyards are 
becoming more and more similar to assembly factories 
using prefabricated elements produced by smaller 
specialised plants. This approach considerably reduces 
costs and, what is more important, the production time, 
which for numerous shipyards may be a decisive factor 
on the present-day competing market. 


Certainly, new solutions bring new problems, which also 
are to be solved for the steel sandwich panels. One of 
basic problems refers to the joints between the panels, 
and between the panels and the classical ship structure. 
Statistically, the most frequently used type of joint is the 
panel-panel heading joint in the longitudinal arrangement. 
From among frequently proposed solutions of this joint 
two types were selected taking into account technological 
realisability, which were the cover plate joint and the 
rectangular profile joint. Parametric models were worked 
out for these two joints. A solution was searched which 
would secure relatively low coefficients of geometric 
concentration at small mass and acceptable deformation. 
The analysis of the obtained results made it possible to 
suggest the optimum joint geometry. 


Steel sandwich panels seem to grow in importance. 
They are likely to take a special place in the production 
of mega-yachts, the demand for which is enormously 
high, as their production making use of plastics is highly 
uneconomical. 
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ABSTRACT 


In this paper anew continuous model for vibration analysis of a beam with an open edge crack is presented. 

A quasi-linear displacement filed is suggested for the beam and the strain and stress fields are calculated. 

The equation of motion of the beam is calculated using the Hamilton principle. The calculated equation of 

motion is solved with a modified weighted residual method and the natural frequencies and mode shapes 

are obtained. The results are compared with those obtained by finite element method and an excellent 

agreement has been observed. The presented model is a simple and accurate method for analysis of the 
cracked beam behavior near or far from the crack tip. 


Keywords: vibration, crack, beam, natural frequency, mode shape, weighted residual 


INTRODUCTION 


Fatigue and cracks are important subjects in industrial 
machineries which can lead to catastrophic failures in certain 
conditions. The importance of the early detection of the cracks 
takes researchers to study various aspects of the behavior of 
a structure defected by cracks. One of these aspects is the 
vibration of cracked structures. Crack development in a system 
changes the vibration behavior. With measurement and analysis 
of these vibrations the cracks can be identified well in advance 
and appropriate actions can be taken to prevent more damage 
to the system. 

The vibration behavior of cracked structures has been 
investigated by many researchers. Dimaragonas presented 
a review on the topic of vibration of cracked structures [1]. 
His review contains vibration of cracked rotors, bars, beams, 
plates, pipes, blades and shells. Two more literature reviews 
are also available on the dynamic behavior of cracked rotors 
by Wauer and Gasch [2, 3]. 

Cracked beam is one of the structural elements which has 
been studied by researchers. There exist three methods for 
vibration modeling of cracked beams: discrete models with 
local flexibility model for crack, continuous models with 
local flexibility model for crack and continuous models with 
continues model for crack. 

For the first time, Dimaragonas suggested the local 
flexibility method for modeling the crack [4]. He assumed the 
crack to be a rotational spring between two healthy parts of the 
beam. The stiffness of this spring obtained from the concept of 
J-integral in fracture mechanics. This local flexibility idea has 
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been followed by several researchers till now. Some researchers 
modeled two healthy half beams discretely and added the 
flexibility of the rotational spring to the flexibility matrix 
of the system [5,6]. While others modeled two healthy half 
beams continuously and used appropriate boundary conditions 
for each part to link them through the rotational spring [7,8]. 
Some other researchers are tried to modify and improve the 
local flexibility model of the crack by adding one or two linear 
springs beside the rotational one [9]. These methods have also 
been extended for beams with more than one crack [10-12]. 
The local flexibility model for the crack is a simple approach 
and has a relative good result in finding fundamental natural 
frequency of a cracked beam. However this method offers 
no solutions for finding the stress at the crack area under the 
dynamic loads, mode shapes in free vibrations and operational 
deformed shape in forced vibrations. 

Another approach to vibration analysis of cracked beams is 
continuous modeling of the crack. Christides and Barr developed 
a continuous theory for vibration of a uniform Euler-Bernoulli 
beam containing one or more pairs of symmetric cracks [13]. 
They suggested some modifications on the familiar stress field 
of anormal Euler-Bernoulli beam in order to consider the crack 
effect. The differential equation of motion and corresponding 
boundary conditions are given as the results. However in their 
model two different and incompatible assumptions have been 
made for displacement and strain fields. Although the accuracy 
of the results in finding the natural frequencies is acceptable 
for some applications their model is not still reliable for more 
accurate analyses such as stress analysis near the crack tip 
under dynamic loading and mode shape analysis. In addition the 


resulted partial differential equation is complex and dependent 
on some constants which are unknown and must be calculated 
by correlating the analytically obtained results with those 
calculated by finite element in each case. Several researchers 
followed the Christides and Barr approach by modifying their 
method and gained some improvements [14-18]. However there 
still exists the inconsistency between strain and displacement 
fields which causes inaccuracy of the results especially in mode 
shapes and stress analysis. 

In this paper a new continuous approach for vibration 
analysis of cracked beams has been presented. The crack is 
assumed to be an open edge crack. A bilinear displacement 
field has been suggested for the cracked beam and the strain 
and stress fields have been calculated. The differential equation 
of motion of the cracked beam has been obtained using the 
Hamilton principle. This partial differential equation has been 
solved with special numerical algorithm based on Galerkin 
projection method. The required constant needed in this 
model can be obtained using fracture mechanics. The results 
of this study are compared with the finite element results for 
verification. 


MAIN IDEA AND ASSUMPTIONS 


The basic assumption in the Euler-Bernoulli bending 
theory for beams is that the plane sections of beam which are 
perpendicular to the neutral axis remain plane and perpendicular 
to the neutral axis after deformation. In the presence of an 
edge crack, the planes will not remain plane after deformation 
particularly at the vicinity of the crack due to a shear stress near 
the crack tip which leads to warping in plane sections. Thus at 
the vicinity of the crack the displacement field is completely 
nonlinear. For the planes far from the crack tip, the warping will 
be smaller and the displacement filed can be assumed linear. In 
order to have a better sense of the bending in a cracked beam, 
a real model has been produced in this research and the mid 
span crack behavior under a pure bending moment can be seen 
in Fig. 1. The beam is made from a linear elastic material with 
low modulus of elasticity and a U-shape notch at the mid-span 
as a crack. 


Deviation line 


Fig. 1. A linear elastic cracked beam subjected to pure bending 


Near the crack area the plane sections will no longer 
remain plane. With a good approximation it can be supposed 
that each plane section turns into two straight planes after 
deformation. The horizontal line passing through the crack tip 
is called “deviation line” in this research which is shown in 
Fig. 1. Each straight plane section turns into two planes with 


different slopes one beneath and the other above the deviation 
line. The slope difference between these two planes decreases 
while the distance from the crack increases. These two straight 
planes connect to each other through a nonlinear part near the 
deviation line. 

In order to find the stress, strain and deformation functions 
for a cracked beam in flexural vibrations a bilinear displacement 
field for the beam has been suggested in this research. In fact 
it is assumed that each plane section turns into two straight 
planes after deformation. The essential assumptions used in 
this research can be listed as follows: 
the beam is slender and prismatic 
the crack is considered to be an open edge notch 
the deformations are supposed to be small 
the plane strain assumption has been used in this research. 
Consequently the displacements along y-axis have been 
neglected 
the stresses are small enough and the crack does not 
grow. 


+ > > 


* 


DISPLACEMENT FIELD 


On the base ofthe above assumptions and explanations, the 
following displacement field is introduced for a cracked beam 
in flexural vibrations: 


w=w(x,t) 
v=0 (1) 
u(x, z, t) =u, (x,t)—zy(x,t)+ A(x, z,HA(z) 


In which u, v, w are the displacement components along 
x, y and z axes. u(x, t) is the longitudinal displacement of the 
deviation line along the x-axis and (x, t) is the slope of the plane 
sections below the deviation line. In an Euler-Bernoulli beam 
theory by neglecting the shear stress effect one has w(x, t) = 
= w (x, t). In a cracked beam the shear stress near the crack 
tip cannot be ignored thus y(x, t) is different from w(x, t). 
However far from the crack the shearing stress decreases 
gradually and (x, t) tends to be equal to w (x, t)h(z) is the 
unit step function which is equal to zero for z < 0 and 1 for 
z > 0. Accordingly the term w (x, t)h(z) can be considered as 
the extra displacement of the plane sections above the deviation 
line. Fig. 2 shows these parameters graphically. 


z A(x,z,) 


Fig. 2. Deformation filed definition of a cracked beam 


The additional displacement of the plane section above 
the deviation line has its maximum value at the crack faces 
and decreases gradually with distance from the crack tip. This 
additional displacement is a nonlinear and complex variable 
with respect to x. Here in this research an exponential regime 
has been assumed for the function A(x, z, t) along the x-axis 
as follows: 
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[x-xel 


A(x, Z,¢) = ọ(z,t)- e 4 sgn(x -x,) (2) 


In equation (2) a is a dimensionless exponential decay rate which will be obtained later in this paper, x, is the crack position, 
d is the depth of the beam and sgn(x — x,) is the sign function which is -1 for x <x, and +1 for x > x,. The application of sign 
function is due to the fact that the additional displacement function has a discontinuity at the position of the crack and the sign 
of its value changes when passing through the crack tip. 

In order to find @(x, t), zero normal stress condition at the crack faces can be used. The normal strain function can be found 
using equation (2): 


[x-xel 


E, EU, =U ZY -Zg te 4 h(z) (3) 


In which the subscript ,x denotes the partial derivative with respect to x. The normal stress at the crack faces should be zero 
so one has: 


d 
251) = — (Wo (os!) — 2 (xa) 4) 


To avoid discontinuity at the crack tip and considering the nonlinearity at the crack tip, the function @(x, t) is modified in 
this paper as follows: 


d - p= 
mz- fu (xX,t)— u (xX, t)e (0) (5) 
pele 0, : 
In equation (5) B is a dimensionless parameter and will be discussed later in this paper. 
EQUATION OF MOTION 


Now the strain field can be extracted from the displacement field. The only nonzero components of the stress field are £, 
and y,, as follows: 


[x-xel 


-pZ “og tel 
E. =ZU =U, -Zy —| u, (x t-u, (x De -zwy (x.t hze 4 
E O,x WY. ast c? O,x E Wy c? 
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1 1 -4 d - 
a a w,-yY+ Ë innte nm Can e 4 sgn(x-x,) 


The normal stress energy of the beam can be obtained using the following relation: 
o1 ol 2 
V =>] oea N Efe dV (7) 
v v 


In which V is the normal strain energy function, V is the volume of the beam and E is the modulus of elasticity. 

In this research the cracked beam is assumed to be slender. So the Euler-Bernoulli assumption can be used and one can 
neglect the shear strain energy in compare with the normal strain energy. Similar to a normal Euler-Bernoulli beam the average 
shear strain in each cross section can be assumed to be zero. So the following relation can be hold: 


[v.44 =0 (8) 
A 
The kinetic energy of the cracked beam can be also calculated as follows: 
1 ; 
T=- | oway 0) 
23, 


In equation (9) similar to the Euler-Bernoulli beam theory the rotational moment of inertia has been neglected. 
Using the Hamilton principle one has: 


t 
Sf (T-V)dt=0 (10) 
to 
Now using equations (6) to (10) and performing appropriate calculations the following equations can be obtained: 
PA 
Ug x = ZY. t KZ, (Xo t)e i 
be (11) 


W y = W x 7 KW o (x, t)e d 
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In equation (11) z is the vertical coordinate of the centroid of 
the cross section and Z, is the vertical coordinate of the centroid 
of the healthy part of the cross section as shown in Fig. 3. The 
parameters k, and k, are geometrical dimensionless constants 
defined in equation (13). 


Fig. 3. A cracked beam parameter definition 


And finally the equation of motion for free vibrations of 
a cracked slender beam can be obtained as follows: 


2 
2 ET, 
ax? 


[x= 


zg ¢ 
[e a TKW Jee * 


}: pAw, =0 (12) 


In equations (11) and (12) the parameters k, , and « are 
geometrical dimensionless constants which can be defined as 
follows: 


k,=—fe add 
A, A, 
A, 
> A +kA 
k, Alf, hk] 
Zh 
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ka= (13) 
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A B, 
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K=k,-— 


(Az Z, (1-k.)-1,, 


Where A is the cross section area of the beam, A, is the 
area of the crack face, I, is the moment of inertia of the crack 
face about y-axis and Lis the moment of inertia of the cross 
section about the horizontal axis passing through the centroid 
of the cross section n. 

The equation (12) is the main result of this investigation. In 
a normal beam the geometrical parameter «K is zero and hence 
equation (12) results in the Euler-Bernoulli beam equation for 
slender beams. The dimensionless exponential decay rates 


(a, B) are the only factors which has not been discussed yet. In 
the next section the parameters a and ß are calculated. 


EXPONENTIAL DECAY RATES a AND B 
CALCULATION 


When a pair of bending moments M are applied to the 


cracked beam an additional relative rotation 0° will exist 
between two ends of the beam due to the crack as shown in 
Fig. 4. 


Fig. 4. Additional remote point rotation of a cracked beam 


It can be shown that for a cracked beam under pure bending 
equation (12) will turn into the following form: 


dw M Ko -al 

L =—_|1-—~e" 4 
a& ET, l+K 

Solving equation (14) will result in the load-deflection 


relation of a cracked beam under static pure bending. The 
results are as follows: 


(14) 


MIx do g A= 
TA Poa =: e E 
a +k 
ee (15) 
M ge Ei K am 
txt, = e x>xX, 
El, 2 a l+K 


In which the constants c,, c,, c, and c, will be as follows: 


c d K oa 
? a? 14K 
Cy =C, +2x a 
d a (16) 
LG ld g am 
C= = 
2 1 l æ 1+kK 
6 =G72— "Eaa 
a |l+K 


Now using equation (15) one can obtain the additional 
remote point rotation 6° as follows: 


g = R (0) = 6(0)) a (O aaa (/ ) z Al )) = 
2M d Kk 


El, a 14K 


In equation (17) the parameter x is a function of B. However 
comparing the finite element results with those obtained by this 
model shows that the parameter ( has a very large value and 
accordingly it can be assumed that the parameter $ is infinity. 
So in equation (17) one can substitute « with limk. 

On the other hand the additional remote point rotation 0° 
has been obtained by empirical methods as follows [19]: 


(17) 
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Equating the right hand sides of cas and (18) will result in parameter 4 values as presented in Fig. 5. 


Exponential decay rate [a] 


0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 
Crack depth ratio [a/d] 


Fig. 5. Exponential decay rate á versus crack depth ratio a/d 


In the next section the partial differential equation (12) has been solved and the natural frequencies and mode shapes have 
been calculated. 


EIGEN SOLUTION 


In order to find the natural frequencies and mode shapes of a cracked beam, the equation of motion presented in equation (12) 
must be solved. However this equation cannot be solved analytically and a numerical method must be used. It can be assumed 
that the solution is a harmonic function so one has: 


w(x,t) = X(x)e™ 12) 
In which o is the natural frequency of the beam. Substituting equation (19) into (12) and assuming EI, to be constant along 
the beam the following eigenvalue problem will be resulted: 


2 gea 
a Xalak 4 eon (20) 


Equation (20) has a special form and contains a singular function and depends on the value of the solution at the crack 
position. Theses anomalies prevent one to use the normal weighted residual solution for this Strum-Liouville problem. In 
anormal Strum-Liouville problem one can easily consider the function X to be in the form of Xc S(x) in which S,(x)are the shape 
functions which satisfy the physical boundary conditions. However in this research the results show that such an approach will 


X-—X6E | 
lead to divergence of the results. Since the functione * ¢ in equation (20) is not a smooth function it seems that the solution 
specially for larger crack depth ratio tends to have large derivatives near the crack tip. Accordingly extracting the value of X"(x,) 
from X by derivation can lead to large fluctuations in the results and divergence. In order to avoid the divergence problem the 
function X” and the value of X"(x,) is not extracted from Xby direct derivation. Instead the X” is discretised independently from 
X and then a constraint equation provided to link X” to X. 
Considering the above discussion the following relations can be written: 
|X-Xe| 


xX" + KX"(x,)e- t= oS (x) 


X = y c; S(x) 
i=l 


In which c, and c’ are two independent set of constants and functions S(x) are the shape functions which must satisfy the 
physical boundary conditions. Substituting equation (21) into (20) and multiplying two sides of the equation by S. (x) then 
integrating along the inet) of the beam one has: 


Ye il S/"(x)S; (x)dx 2 2 Ye ai S/(x)S,(x)dx=0 jf =1,2,...,N (22) 


EI, 


(21) 
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Or in the matrix form: 


pA» 
K; |la]--@ LF, 


Tel =0,K,= f S (5S, (x)dx, P, = [ S,(x)S, (x)dx 


(23) 


On the other hand from equations (23) the following relation can be obtained: 


N , n 1 A 
ay =— ‘Oy 
es O= Yes 


(24) 


Multiplying two sides of equation (24) by S(x), integrating along the length of the beam and writing the equations in the 


matrix form one has: 


-l 


la]; f S (x)S,(x)dx R, -— f 50S, (de 


[a]=[2;] [7 25) 
Substituting equation (25) into (23) the following equation will be resulted: 
på , 7 = 2 al 
K; Ea [M, | [e,]= 0, [M, | = P, | Q; R; (26) 
, 


The natural frequencies and corresponding mode shapes for the cracked beam can be calculated solving the matrix eigenvalue 
problem of equation (26). In the next section the results are presented for a simply supported beam with rectangular cross- 


section. 


RESULTS FOR A SIMPLY SUPPORTED BEAM WITH RECTANGULAR CROSS-SECTION 


In this section the new approach has been applied for free vibration analysis of a simply supported slender prismatic cracked 
beam with rectangular cross-section. In such a beam the exponential decay rate B can be assumed to be infinity and the exponential 
decay rate a can be calculated from equations (17) and (18). The geometrical factor x and the exponential decay rate a are as 


follows: 


In a simply supported cracked beam the shape functions 
S,(x) can be assumed to be in the form of sin(iđx/1) which satisfy 
the physical boundary conditions. The natural frequency and 
mode shapes can be calculated using eigenvalue problem of 
(26). In this research the number of shape functions N is set to 
be 100. In order to generalize the results the natural frequencies 
of the cracked beam have been divided to the corresponding 
values for a normal beam. Figures (6), (7) and (8) show the 
fundamental, second and third natural frequency ratios of the 
cracked beam respectively. In figures (6) to (8) the natural 
frequency ratios have been plotted versus the crack depth ratio 
(a/d) for several crack positions. 

In figures (6) to (8) the results of finite element (FE) 
analysis are also presented for verification. The finite element 
results have been obtained using ANSYS software. In order to 
have an accurate and reliable model the PLANE183 singular 
element has been used in the cracked area [20]. This element 
is an 8-node quadratic solid singular element which specially 
designed for crack analysis. In this research a fine mesh has 
been used at the vicinity of the crack and dependency of the 
results to the mesh size has been checked. In all of the results 
there is a good agreement between analytical results and those 
obtained by FE analysis. 

As it can be seen in figure (6) the reduction rate of the 
fundamental natural frequency has a direct relation with the 
position of the crack. This rate reduces for cracks which have 
more distance from the mid span of the beam. For the cracks 
at x/l = 0.1 the fundamental natural frequency only drops 
nearly | percent when the crack reaches to the half of the beam 


SSS SSS SSS SS eee 
593| © |~-1969| ©] +3714 ©] -3584/ E +1312] 2 
d d d d d 


off 


(27) 


depth while for the cracks at the mid span this value is about 
15 percent. 

The dependency of the reduction of the natural frequency 
to the crack position is also seen in the first few natural 
frequencies. For the cracks at the mid span the second natural 
frequency remains nearly constant with the crack depth because 
this point coincides with the node of the second vibration mode 
of the beam. 

Figures (9) to (11) show the first three normalized mode 
shapes for a cracked beam with a/d = 0.5 and x /l = 0.1, 0.3, 
0.5. Comparison of the analytic and finite element results in 
this set of figures shows the efficiency of the model presented 
in this research. 
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Fig. 6. Fundamental natural frequency ratio of a cracked beam versus 
crack depth ratio (a/d) 


POLISH MARITIME RESEARCH, No 2/2008 37 


= 
$ — Analytical results for x /=0.1 
T —— Analytical results for x A=0.2 
3 — Analytical results for x =0.3 
Bo» —— Analytical results for x /=0.5 
¢ + FE results for x A=0.1 
IE os a FE results for x A=0.2 
= FE results for x =0.3 
os © FE results for x A=0.5 
ka) 04 02 03 os 09 1 


os os os 
Crack depth ratio [a/d] 
Fig. 7. Second natural frequency ratio of a cracked beam versus crack 
depth ratio (a/d) 
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Fig. 9. First three normalized mode shapes of a cracked beam with x /1=0.1 


and a/d=0.5 (- ): Analytical results; (¢*¢*): Finite element results 
= 04 
3 03 
8 02 
Bot 
Zo 
Bos 
G02 
E 03 
Za 01 02 03 04 os os o7 o8 os 1 


Position ratio [x/l] 


Fig. 10. First three normalized mode shapes of a cracked beam with x / 
1=0.3 and a/d=0.5 (- }: Analytical results; (¢*¢*): Finite element results 
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Fig. 11. First three normalized mode shapes of a cracked beam with 
x/=0.5 and a/d=0.5 (—_): Analytical results; (****): Finite element results 
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CONCLUSIONS 


O Anew continuous theory for flexural vibration analysis of 
a beam with an edge crack has been introduced in this paper. 
The crack is assumed to be an open edge crack. A bilinear 
displacement filed has been suggested for the beam and the 
strain and stress fields have been calculated. The extraction 
of the strain field from the displacement field is based 
on elasticity rules and hence the displacement and strain 
fields are completely compatible. The governing equation 
of motion has been obtained using the Hamilton principle. 
The required constants in this model are calculated from 
empirical formula in fracture mechanics. The equation of 
motion has been solved for natural frequencies and mode 
shapes using a special numerical algorithm presented in 
this paper. 


O The analytical results have been compared with finite 
element results and an excellent agreement has been 
observed. The model is accurate for both natural frequencies 
and mode shapes calculations. Results show that for 
a simply supported beam the reduction of natural frequency 
is a function of crack depth ratio as well as the crack 
position. 


O The presented model is a simple and accurate model which 
predicts the behavior of the cracked beam and its results are 
reliable near the crack tip and far from it. This model can 
be used for dynamic stress and strain calculations near the 
crack tip. 
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Determining strength of the modified wood 
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ABSTRACT 


The article presents results of investigations of the effect of pinewood modification with methyl 
polymethacrylate on the anisotropy of its strength in complex stress conditions. The investigations aimed 
at determining the resistance of the examined wood to stretching, compression, stretching with shear, and 
compression with shear. For the investigations oriented on shear, shear with stretching, and shear with 
compression, a special specimen was prepared which differs by notch geometry from a typical Iosipescu 
specimen. A new test machine is described in the article, which is equipped with special specimen holders 
to perform investigations in complex stress conditions. Crack patterns recorded for the natural and modified 
wood are presented. For all tests, numerical FEM simulations were performed to obtain stress distributions 
inside the specimens. The calculated stress distributions were visualised as contour line projections for 
the natural and modified wood. 


Keywords: wood modification, wood-polymer composite, Iosipescu specimen, 
stress distribution in the specimen under load. 


INTRODUCTION 


Numerous authors [5, 6, 10, 11] claim that the properties 
of the modified wood relating to its strength and deformation 
increase or decrease in proportion to the amount of the 
modifying substance. 

All changes of mechanical properties of the wood exposed 
to the modification cannot be fully predicted without relevant 
experimental tests. 

The goal of the reported activities was to determine the 
effect of pine whitewood modification of its strength anisotropy 
in complex stress conditions. 

Wood is an orthotropic material revealing remarkably 
different strength characteristics depending on the anatomic 
direction. The literature makes distinction between three wood 
directions: longitudinal, tangential, and radial. The adopted 


Longitudinal 


Annual ring 


Tangential 
seo) 


Radial R, x; 


Fig. 1. Orientation of anatomic directions 
of the wood in rectangular coordinate system 
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coordinate system is shown in Fig. 1 — the L axis is in the 
anatomic direction along the fibres, the R axis is in the radial 
direction, and the T axis is directed tangentially to the surfaces 
of particular fibre layers. 


COURSE OF INVESTIGATIONS 


In order to examine the effort of the examined materials, 
tests were performed to determine wood resistance to stretching, 
compression, stretching with shear, and compression with shear. 
The shapes and dimensions of the specimens used for tests were 
taken from relevant standards [22+25]. 

The specimens were formed from a wooden square timber, 
seasoned in natural conditions in the temperature of 295+2K. 
The average humidity of the specimens was equal to 8%. 
The tests were performed on natural wood specimens and 
those exposed to the modification process. The modification 
consisted in saturating the wood with the hydroquinone methyl 
ether modified methacrylate (MM), with further thermal 
polymerisation. The degree of specimen saturation with the 
monomer was controlled by adjusting the time during which 
the specimen was in the autoclave, in the way described in 
[21]. The modified wood with different contents of methyl 
polymethacrylate was marked, respectively, as K0.0, K0.35, 
K0.43, K0.48 and K0.56, where the numbers represent the 
number of kilograms of methyl polymethacrylate with respect 
to 1 kg of dry wood. 

The stretching and compression tests were done on 
a universal test machine MTS-8 10.12, while the shear tests were 
performed on a machine designed and produced especially for 


this purpose [21]. Wood resistance to stretching, compression, 
and shear was determined based on the arithmetic mean of 6 
to 10 specimens. The mean values and the statistic processing 
are discussed in [21]. The results of the mean value analysis 
are given in Table 1. 


Tab. 1. Resistance limits to stretching, compression, 
and shear of natural and modified wood 


Type of 
material 


where: Ra. - resistance to stretching along fibres; R, - resistance to 
compression along fibres; Ras - resistance to stretching across fibres; 
Ro - resistance to compression across fibres; R, - shear resistance 


RESISTANCE OF THE MODIFIED WOOD 
TO IN-PLANE SHEAR 


Prognosing strength of such composite layer structures 
as the wood is of vital importance in designs making use of 
this material. An important parameter which should be known 
in wood-based designs is wood resistance to in-plane shear. 
The literature sometimes gives contradictory data on the 
levels of wood resistance, which is not surprising [15]. The 
absence of reliable methods measuring material’s resistance 
to in-plane shear results in little confidence in official design 
limitations, a consequence of which is taking into account 
excessive margins in design activities involving in-plane 
shear strength. 

Methods of composite resistance investigations can 
be divided into tests making use of tubular specimens and 
plane specimens. Preparing wooden tubular specimens for 
investigations in complex load conditions is difficult due to 
the anatomic structure of the wood (annual rings), therefore in 
the reported case a decision was made to perform tests using 
plane specimens. 

So far, a specimen which was found to be most useful, due 
to its shape, for composite tests is the Iosipescu specimen [4, 
8, 9]. Attempts to adapt it to testing composites in plane stress 
conditions, following the method based on the Arcan instrument 
[4], were initiated by Broughton [7]. Its applicability for shear 
tests of various layer composites was confirmed by numerous 
investigation reports [1, 2, 13, 17=19]. 

The failure of the Iosipescu specimen can be reached in the 
uniform stress conditions, although the presence of undesired 
transverse stress makes the interpretation of the shear resistance 
more complicated. 

Below, assumptions are presented which refer to the 
measurement of the shear resistance of the composite Iosipescu 
specimens in the test bearing the name of the V method (ASTM 
D5379M:93). The definition of the shear resistance is, generally, 
not clear [14, 21]. For instance, for the same material Broughton 
[7] says about the “failure stress” equal to 58 MPa, while 
Morton [12] says about the “shear resistance” equal to 68 MPa. 
At the same time Adams [1] says about the “shear resistance” 
equal to 115 MPa for the same material and specimen, but 
another definition of the failure point [14, 15, 21]. 

When testing plane specimens we have to take into account 
that reaching the failure criterion generates the necessity for the 
appearance of shear. An important feature of the test is stress 
uniformity in the examined specimen area. 


Fig.2.: a) modified losipescu specimen, b) specimen with summer wood 
upper layers, c) specimen with spring wood upper layers. p = 45°, g =5 
mm - specimen thickness, b = 20 mm — specimen width, p = 5 mm — length 
of measuring section, L = 80 mm — specimen length 


Preliminary tests performed ona typical Iosipescu specimen 
(with an alternate system of layers, not always parallel to each 
other, sometimes twisted) have revealed that correct results 
cannot be obtained. During shear tests of this piece, it is not 
only the shearing stresses which are observed at the bottom 
of the notch. When the specimens are loaded with transverse 
load, their bending is observed, with the resultant normal 
stress generated by the bending moment. The wooden fibres 
are deformed (bent and stretched). Recognising the beginning 
of specimen failure is difficult [21]. 

Taking into account the experience gained by Iosipescu [8], 
a specimen was prepared which was then tested many times 
to reach the final shape and dimensions shown in Fig. 2. This 
specimen differs considerably from the typical Iosipescu piece 
by the notch geometry. Adopting such a specimen geometry is 
aresult of long lasting experiments, numerical calculations and 
observations recorded during the tests. The depth of the notch 
with respect to the thickness of layers was selected in such 
a way that the beginning of the notch is at its edge [21]. 

The examined specimens had various configurations of 
softwood and hardwood layers. When preparing the specimens, 
care was taken that the wood layers were situated parallel to 
each other and had constant thickness along the specimen 
length. 

There are test machines, produced in series, which are 
equipped with special holders for performing investigations 
in complex load conditions. However, due to high costs of 
these machines and specific nature of investigations, similar 
machines are frequently designed and produced individually in 
laboratories in which the investigations are preformed. 


F- Pressing force generated 
by the test machine 


Fig. 3. Machine for generating plane stress distribution: 

1 - casing; 2 - console; 3 - base; 4 - guide; 5 - pilot sleeve; 6 - right holder; 
7 - left holder; 8 - insert for left holder jaw; 9 - insert for right holder jaw; 
10 - left locking block; 11 - right locking block; 12 - guide sleeve; 

13 - 912 ball and ball fixing; 14 - M6 bolt. 


The designed and produced machine [21], shown in Fig. 3, 
has made it possible to generate the plane stress distribution, 
Fig. 2. Compared to a traditional test machine [7, 20], the 
position of the specimen axis with respect to the direction of 
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force F can be changed. As a result, depending on the specimen 
position angle with respect to the load direction, stretching or 
compression can be executed along with shear in the central 
part of the specimen. 


Moving / 
blocks \ 


Fig. 4. Scheme of kinetic load of the specimen 


The moving and fixed blocks are preliminary clamped in 
such a way that the Fx, component of force F is smaller than the 
friction force between the blocks and the specimen. However, 
this clamping generates a preliminary stress. Taking this into 
account, the clamping force was selected experimentally at 
a level securing fulfilment of its function and, at the same 
time, not affecting considerably specimen’s strength. Then 
the moving blocks were moved in the direction inclined at 
the angle 6 to the specimen axis. During the tests the value of 
force F was measured in relation to specimen displacement. 
When 6 = 0, the stress is close to pure shear. The test machine 
makes it possible to change the angle within the range -45+45°. 
Specimen stretching and compression in the direction of the 
X, axis was done when the specimen was mounted directly in 
the machine jaws. 

The scheme shown in Fig. 4 reveals some displacements 
of the upper blocks with respect to their lower counterparts. 
The scale of this displacement was set experimentally. It was 
large enough to generate united displacements at the contacting 
edges between the blocks and the specimen, which was of high 
significance for the operation of the FEM model used for the 
numerical stress analysis. 

For 6 = 0, shear of the specimen is observed as provoked 
by the applied vertical displacement U, , Fig. 4. The 
machine that loads the specimen records, ‘via a relevant 
measuring instrument, the force along with the corresponding 
displacement. During specimen shear, the displacement U, 
of the machine right-hand side, corresponding to force F = F, 
is directed along the x, axis. For negative (counterclockwise) 
specimen rotation angles in the test machine, see Fig. 3, force 
F generates shear F, and compression F, within the specimen, 
while for positive (clockwise) angles force F generates shear 
F and stretching F . 

Fig. 5 shows a sample relation between the loads of the 
specimens made of natural and modified wood, and the holder 
displacements for ò= 0°. Numbers 1+4 denote successive crack 
stages of the specimen made of natural wood. The specimen 
failure does not take place immediately after the appearance 
of the first crack, as the arrangement of the wood layers with 
respect to the direction of force action changes at that moment. 
Further increase of load generates next cracks. 

Some authors suggest that the final crack takes place, for 
plastics, in the uniform shear conditions. Unfortunately, no 
investigations were performed to verify this thesis [3, 7, 16]. 
Other authors propose to interpret the shear stress recorded at 
the appearance of final crack as the shear resistance. 

Here, for both the natural and modified wood the shear 
resistance was assumed as that observed at the first crack of 
the specimen, Fig. 5 [21]. 
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Fig. 5. Load - displacement curve with characteristic points of crack 
beginning (1, 2, 3, 4 — successive stages of specimen crack) 


Fig. 6 shows crack patterns observed for the natural and 
modified wood. When exposed to load, the layers of the 
natural wood resembled significantly deformed “threads”, 
and recording the time when the last layer is damaged was 


Fig. 6. Specimen crack patterns: a) natural wood, b) composite 
D-PMM K0.35, c) composite D-PMM K0.56 
(1, 2, 3, 4 — successive stages of wooden specimen crack) 


difficult. Specimens of the modified wood revealed higher crack 
resistance. When the content of the polymer in the composite 
was higher, the material became more resistant to failure, Figs. 
5 and 6c. The nature of crack of the modified wood K0.35 
resembled an intermediate crack pattern between that observed 
for brittle and ductile material. 

The above course of specimen crack referred to the case 
of shear. In a similar way, depending on the selected angle 
of specimen axis with respect to the direction of load, forces 
F were recorded assuming that the specimen resistance 
corresponds to the first crack (point | in Fig. 5). By changing 
the specimen position angle by 15° within the range -45+45° 
the relation was obtained between the force F corresponding 
to the first crack and the angle ò. 

Fig. 7 shows loads of specimens made of natural and 
modified wood in relation to the angle 6 until the appearance 
of the first crack. 


1300 KO0.0 
KO0.35 
FIN] K0.43 
KO0.48 
1100 KO0.56 


900 S 


6 [°] 
60 
-45 -30 -15 0 15 30 45 
Fig. 7. Dependence of specimen load on the angle 6 until the appearance 


of the first crack 6 — specimen rotation angle with respect 
to the horizontal position 


ANALYSING TEST RESULTS 


For all tests, FEM-based numerical simulations were 
performed to obtain stress distributions inside the pieces. In 
the calculation model the plane stress was assumed, and the 
discretisation was done in such a way that each layer of the soft 
and hard wood was modelled in the perpendicular direction by 
two nine-node elements. Physical properties of particular layers 
were described using a generalised Hooke law for orthotropic 
materials. Based on the obtained results, the flexibility matrices 
were determined for the soft and hard layers of the natural wood 
K0.0 and the modified wood K0.56 [21]. 

On average, the specimen consisted of seven softwood 
layers and eight hardwood layers, or the opposite. The ratio 
between the thicknesses of hardwood layers to softwood 
layers was approximately equal to 0.5. The depth of the notch 
with respect to the layers was selected in such a way that the 
beginning of the notch was at the edge of the layer and not in 
the middle of it. The layers were assumed to be situated parallel 
to each other and have constant thickness along specimen 
length. In the measuring section the top layer can be the either 
softwood or hardwood layer, Figs. 2b, c). 

The specimen load was executed by moving the edges 
which were in contact with the moving blocks in the direction 
inclined by ô to the specimen axis. When the force calculated 
from stresses for section x, = 0 was equal to that measured, the 
obtained stress conditions were assumed as corresponding to 
specimen failure. Comparing edge displacements assumed in 
the calculations with the measured values made the basis for 
verification of the calculation model. 


a) 


o[MPa] 


10.02 
O: 


6.90 


o[MPa] 


66 | y=0 033 


-2.46 


750 8.50 9.50 10.50 11.50 12.50 


xı [mm] 


Fig. 8. Distribution of stresses in the measuring sections of the specimens 
made of: a) natural wood K0.0, b) modified wood K0.56 exposed to complex 
stress conditions for ò = 0° 


Fig. 8 presents stress distributions in the measuring sections 
of the specimens made of natural wood K0.0 and modified 
wood K0.56 for 6 = 0. Here, the shearing stresses o,, take 
the highest values, while o,, and o, = 0.250,, are close to 
zero. The stress distributions along the examined measuring 
sections of the specimens made of natural and modified wood 
differ between each other. For the natural wood, the saw-like 
shape of the diagram suggests the appearance of increased 
stresses in the hardwood with respect to those observed in the 
softwood. In the modified wood, on the other hand, the stress 
distribution is close to uniform, as elastic properties of the 
layers are similar to each other. Modifying the soft layers led 
to the “homogenisation” of the material, as a result of which 
the stress distribution curve is more uniform. 

Figs 9 and 10 show, in the form of contour lines, the 
distributions of stress 6,, calculated using the FEM method 
for 56 = 0 for natural and modified wood, respectively. The 
distribution of the tangential stress in the measuring section of 
the specimen made of natural wood, Fig. 9, reveals much higher 
effort of the hardwood than that of the softwood. Moreover, 
within the entire measuring section the material effort is more 
than twice as high as in the corresponding section of the 
modified wood, see Fig. 10. This effect is believed to have 
been caused by the “homogenisation” of the material in the 
hardwood layer direction [21]. 

For the specimen made of natural wood, stress distributions 
in the analysed section are shown Figs. Ila and b, while for 
the modified wood - in Figs. 11c and d. As can be noticed, 
the stresses for the two angles differ not only by a sign, but 
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also by an absolute value. These differences are much higher 
for the natural wood than for the modified wood. Similar 
differences in stress distributions in the measuring sections 
of the specimens made of wood K0.0 and K0.56 were 
recorded at complex load conditions for 6 = 90 and 6 = -90. 
In the specimens exposed to load for 6 = 90 and 6 = -90° the 
dominating stresses o,, are accompanied by limited stresses 
o, and negligible shear stresses. For the both cases (6 = 90 


Fig. 9. Distribution of o,, stress in selected natural wood specimen area 
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and ò = -90°) the tests were performed on a universal test 
machine MTS 810.12. 

The shear resistances and stress distributions were determined 
in complex load conditions (shear with stretching, and shear with 
compression) in the measuring sections of the specimens made 
of natural and modified wood. If the top layer in the measuring 
section was the hardwood layer, the stresses were higher than 
in case when the softwood layer was the top layer. 


X] 


Fig. 10. Distribution of o, stress in selected modified wood specimen area 
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Fig. 11. Distribution of stresses in the measuring sections of specimens made of natural wood K0.0 at complex stress conditions for: 
a) 6 = 45°, b) 6 = -45° and modified wood K0.56 for: c) 6 = 45°, d) 6 = -45° 
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FINAL CONCLUSIONS 


The performed tests have made it possible to formulate the 
following conclusions: 

O When being compressed along the fibres, the effort of the 
modified wood can be nearly twice as high as for the natural 
wood. In case of compression across the fibres this ratio 
is even higher and can exceed 4. Also during stretching, 
remarkable increase of strength of the modified wood is 
observed. 


O Stress distributions in specimens made of natural and 
modified wood and exposed to load reveal differences in 
levels and shapes. These differences are much higher for 
natural wood. 


O When exposed to load, layers of the natural wood resembled 
significantly deformed “threads” and capturing the time 
when the last layer was damaged was difficult. The 
specimens made of modified wood revealed higher crack 
resistance. Increasing the content of polymer resulted in 
“homogenisation” of the wood structure, as a result of which 
the material became more resistant to failure. 


O Inthe elements made of natural wood remarkable differences 
in the efforts of soft and layers were observed, while in the 
superficially modified wood the efforts of the soft and hard 
layers situated close to the surface were almost the same. 


O In case of bent or twisted elements, a sufficient approach 
is to modify the wood only superficially. The softwood 
layers are strengthened and take over some part of the load 
from the hardwood layers. This way the stress distribution 
becomes more uniform. 
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ABSTRACT 


This paper devoted to ship combustion engine diagnostics, presents short characteristics 
of experimental test techniques applicable to state assessment of piston and gas turbine 
engines presently being in use in Polish Navy. Some metrological and organizational aspects 
of realization of necessary diagnostic tests as well as measurements of control parameters 
of engines of low susceptibility to control, were highlighted. Such engines are installed 
a.o. on board American origin frigates of Oliver Hazard Perry type as well as Norwegian 
origin submarines of Kobben type which have been recently incorporated into Polish Navy 


(i.e.General Electric turbine engines, Detroit Diesel and Mercedes Benz piston engines), as well as older 
engines of Russian production (Zorya turbine engines, Zvezda piston engines). Applicability of particular 
diagnostic methods was also assessed from the point of view of their effectiveness in identifying engine 
failures in conditions of its operation in ship power plant. Recently realized tests have made it possible to 
conduct operation of all piston and turbine engines used on board Polish Navy ships in accordance with their 
current technical state. On the average 120-130 diagnostic expertises are yearly issued, and the carried-out 
tests cover main functional systems of engines together with their control and safety subsystems. 


Keywords: marine diesel and gas turbine engines, technical diagnostics, base diagnostic system, operational examinations 


INTRODUCTION 


Main research activity of the Institute of Ship Construction 
and Operation (IKiEO), Polish Naval Academy (AMW), is 
focused on the problems of diagnostics of ship combustion 
engines (both of piston and turbine kind). The activity has 
been developed continuously since 1982 owing to efforts of 
the following research teams: that dealing with turbine engines 
and supervised by Prof. Charchalis, A., that dealing with 
piston engines under supervision of Prof. Piaseczny,L., and 
Dr Polanowski, S., as well as that dealing with diagnostics of 
flow systems of turbine and piston engines under supervision 
of this author. 

Introduction to Polish Navy, in recent years, driving 
engines of a new type of low susceptibility to control, makes 
it necessary to search for new, alternative diagnostic methods 
providing possibility of carrying out a complex assessment 
of technical state of the engines independently from very 
expensive servicing. 

In this paper is described the range of scientific research 
undertaken in the Institute, which have been traditionally 
focused on the area of operation of ship machines and devices, 
including a.o. development of assessment methods of technical 
state of ship engines (and not only) being in service conditions 
on board a ship. Also, selected scientific problems dealing with 
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diagnostics of ship piston and turbine combustion engines, 

being recently under research in the Institute, are presented. 
The below presented illustrative material was divided in two 

independent parts: the first -devoted to diagnostics of piston 

engines, the second— to diagnostics of turbine engines. The 

described “cuisine” of diagnostic techniques contains selected 

methods of diagnostic activity based on the following scope 

of measurements of: 

Gasodynamic parameters of working medium, of high- and 

low variability, 

Vibrations and their spectral and correlation analysis, 

Metallic contamination in lubricating oil, 

Torque, 

Emission of toxic compounds contained in exhaust gases, 

Parameters of supplied fuel, 

As well as endoscoping examinations. 


HHHHHH 


It should be added that a part of diagnostic instruments 
can be sucessfully applied to assessment of technical state 
both piston and turbine engines, e.g. an endoscopic set or gas 
analyzers. However procedures for realization of diagnostic 
tests are different as they are adjusted to a given type of engine, 
from case to case. A strategic aim of the presented diagnostic 
methods, continuously developed and modernized, is to bring 
them to practice of operation of ship engines installed on Polish 


Navy ships and subjected to diagnostic survey in accordance 
with their current technical state. 


DIAGNOSTICS OF PISTON ENGINES — 
A HISTORICAL OUTLINE 


Since 1982 in the Institute a series of scientific research 
investigations concerning the diagnosing of piston combustion 
engines has been realized. The undertaken scope of work 
has inspired actions leading to elaboration of effective and 
complex diagnostic methods applicable in service conditions 
of ship engine. A strategic aim was to pass to engine operation 
process in accordance with its technical state. As a result the “ 
Basic diagnostic system” was implemented into the operational 
system of Polish Navy ships. Its concept and realization resulted 
from: 

# the task project no. 148-39/C-SO/93.: „Basic system for 
diagnosing and predicting the technical state of piston 
combustion engines”, 

# the research project no. 9T12D00211.: ,, A method of 
diagnosing the main design systems of piston engines 
and compressors with application of analysis of vibration 
envelopes in the domain of crankshaft rotation angle”, 

realized by the team under supervision of Dr Polanowski. 

(Fig. 1). 

Since 1992 have been systematically carried out diagnostic 
tests of piston engines of every type appplied to propelling 
Polish Navy ships. Periodic assessment of technical state of the 
engines is realized by an experienced research team specialized 
in application of the newest methods of diagnosing in difficult 
conditions of ship operation. Realization of such tests is 
based on successive indication of engine cylinders during 
its steady operation under load within representative ranges. 


NIP max dĦIP = 0.3 POOH max APCON = 1. 
Kbar] nin aHIP = 9.3 [bar] nin APCO = 3. 
15.2 116.3 


12.5 deg 
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With the use of special analyzers of fast-varying quantities 
[20] is carried out comparative, statistical and subject-matter 
analyses of developed indicator diagrams, runs of vibration 
accelerations generated by working mechanisms connected 
with cylinder system and transmitted to a measurement point on 
engine head, as well as runs of other quantities characterizing 
working processes in engine cylinders, recorded in the domain 
of crankshaft rotation angle. The so determined diagnostic 
measures: mean indicated pressure, indicated power, maximum 
combustion pressure and rate of cylinder pressure increase dp/ 
da etc provide important information on overall technical state 
of constructional elements of combustion chamber. Results 
of the tests are used for analyzing trends, assessing changes 
in technical state of engines and preparing an operational 
decision. 

The presently performed investigations make it possible 
to carry out, in accordance with the technical- state control 
procedure, operation of about 100 engines installed on board 
Polish Navy ships. On average, 75—80 expertises on technical 
state of the engines are yearly issued [21,23]. And, only those 
engines are subjected to the diagnostic control system which 
are fitted with indicator valves. 


ENDOSCOPIC DIAGNOSTICS 


Novel ship engines are fitted with more and more perfect 
control — measurement systems to measure parameters 
characterizing their loading state. Despite that, from service 
practice are known cases of severe failures of engines, whose 
primary causes have not been found in a proper time, e.g. 
failures due to excessive vibration resulting from loss of 
stability of mechanical system and - in consequence - occurence 
of a resonance. 


Fig. 1. The diagnostic team of Dr Polanowski at the stand for testing the cylinder systems of Henschel engines — cylinder pressure analyzer, 
example results of testing the engine combustion chamber 
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Difficulties in recognizing failures of ship combustion 
engines on the basis of thermal-gas-dynamic parameters 
characterizing working processes, are associated with proper 
interpretation of symptoms of a given failure. The symptoms 
are often identified as those arising from natural ageing process, 
contamination or wear of engine elements, determined by 
duration time of engine operation. In such cases external 
symptoms are usually concurrent and difficult for unambiguous 
determination. A special case forms the problem of analysis 
of diagnostic parameters for assessment of technical state 
(tightness) of the PRC system (Piston-Ring-Cylinder) of ship 
multi-cylinder piston engine deprived from possible application 
of cylinder indication. It often happens that trend curve of 
exhaust gas temperature changes only slightly deflects whereas 
piston head becomes burned through due to an injector failure. 
In such cases the only possibility to obtain any diagnosis 
of technical state of engine is ensured by visual inspection 
of its internal spaces with the use of endoscopes (elastic 
fibrescopes and stiff borescopes). In this non-invasive, very 
fast, inexpensive and unambiguous way all doubts bothering 
engine operator can be removed. 

To operation process of ship combustion engines new 
diagnostic test methods are commonly introduced. Endoscopy 
which has been earlier used mainly for assessing technical state 
of flow part of turbine engines is dynamically developing and 
now it serves as a very useful and even indispensable tool for the 
diagnostic team of the Institute in its activity concerning piston 
engines especially those being in failure states or in the case 
of necessary prolongation of time interval between overhauls 
[23]. The team in question is equipped with the endoscopic set 
of OLYMPUS and STORZ firms (Fig.2). 


Set of borescopes 
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Valves in engine head 
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The set is composed of three borescopes differring to 
each other with length of optic part, diameter and observation 
angle of diagnosed elements, namely: 90cm/8mm/90°, 
30cm/4mm/0°, 30cm/10mm/120°. The latter one is especially 
useful in diagnosing combustion chambers and valve seats. 
The endoscopes show their high usefulness in hard accesible 
places of engine, e.g. combustion chambers - especially in 
the case when the disassembling of the head is difficult and 
time-consuming, in turbocharging system or internal spaces 
of mechanisms coupled with engine’s crankshaft. In Fig. 2. 
a way of introducing the endoscope tip to internal space of 
piston engine, is presented. After disassembling the injector the 
endoscope makes it possible to assess technical state of piston 
head, cylinder liner surface and other subsystems mounted 
in it such as: spray nozzles of remaining injectors, inlet and 
outlet valves, start-up valves etc. The endoscopic inspection 
method plays special role in diagnosing the multi-block and 
multi-cylinder engines. For instance, in the engines of radial 
arrangement, of M503 type (having 42 cylinders) or M520 
one (having 56 cylinders), installed in engine room, access 
to lower monoblocks and lower parts of reduction-reversible 
gear is very difficult. In the case of necessary overhauls, such 
engine together with the gear must be coupled out of propeller 
shaft, then heeled and lifted and sometimes even rotated inside 
the engine room so as to ensure access to the 1* or 7 cylinder 
block. From service practice it results that a fibresoscope of 
a sufficiently long light pipe of optical system makes it possible 
to avoid such troubles and thus to save realization time of 
overhauls and their costs even by 25-30% [23]. 

Systematic endoscopic tests carried out in the frame of 
periodical maintenance surveys of ship engines operating in 


Results of endoscopic examinations 


Opened inlet valve 


ra Rais“ 


Erosion pit on piston head 


Fig. 2. Endoscopic examination of a piston engine — endoscopic set, results 
of examination of ship engines, access to internal spaces of a cylinder liner 


Polish Navy have confirmed the high effectivenes of the method 
at relatively easy use of the testing instrumentation. As a result 
of such overhauls were disclosed many defects which — in the 
case of their further building-up - could highly endanger engine 
reliability. Selected defects of ship engines identified during 
endoscopic tests in service conditions have been systematically 
published in scientfic journals, a.o. in [8,14,23]. 


ALTERNATIVE METHODS FOR 
DIAGNOSING PISTON ENGINES 


The scientific research work carried out in the Institute (under 
supervision of this author) is aimed at thorough modernization 
and further development of the so far functioning and proved 
diagnostic system by adding to it elements for controlling the 
main functional systems of engine (i.e. the piston-cylinder 
system, turbocharging system, fuel supply system as well as 
bearing system) in the situation when to assess technical state 
of engines not fitted with standard indicator valves becomes 
necessary. Just such engines (M401, Detroit Diesel and 
Mercedes Benz engines) are installed a.o. on the Polish Navy 
base trawlers of F207M type, American-origin frigates of Oliver 
Hazard Perry type and Norwegian-origin submarines of Kobben 
type, which have been lately incorporated to Polish Navy. 

In the frame of realization of the following research 
projects: 
> „Influence of adjustment of fuel supply system of ship 

piston combustion engine on spectrum of shaftline torsional 

vibrations ”, 

> „Method of using measurements of cylinder internal 
pressure during idle operation to assessment of technical 
state of bearing system of ship piston combustion engine 
in service conditions”, 

> „Diagnostic tests of supercharging system in the aspect 
of changes in combustion process occurring inside engine 
cylinder ”, 

> „Research on influence of operational factors on 
performance and efficiency of propulsion systems of Polish 

Navy ships”, 
> „Technical state assessment of ship engines of the series 

of AL, ASV and M520 types on the basis of measurement 

data collected by main propulsion monitoring systems of 
the ships of 206FM and 660 type” 


was elaborated a number of new diagnostic stands which make 
it possible to indirectly assess the fuel combustion process in 
engine cylinders when their indication is not possible. 


DIAGNOSTICS OF TURBOCHARGING 
SYSTEMS 


Turbocompressor constitutes a sensitive and — according to 
statistics - very unreliable element of four-stroke engine fitted 
with pulsatory supercharging system. Such situation results 
from specific operation conditions of ship engines, which are 
characterized by work at low and variable loads. It leads to an 
incomplete, and not quite thorough, fuel combustion in cylinders 
and further gas-dynamical and mechanical consequences of the 
phenomenon in the form of erosion, decreasing performance 
and lowering efficiency of turbocompressor, loss of stability 
of mechanical system, and at last to vibrations which generate 
accelerated wear of bearings and fatigue cracks of rotor’s 
blades. 

In diagnostic tests of pulsatory supercharging systems 
of ship engines it becomes necessary to identify run of 
changeability of being —at- disposal energy of pressure pulses 


of exhaust gases delivered to turbocompressor during one 
working cycle. On the basis of measurements of instantaneous 
values of the ram pressure p1* and instantaneous values of the 
static pressure p1, p2, recorded at two control cuts, 1 and 2, of 
exhaust gas outlet channel, located at the distance L to each 
other — Fig.3, one determines the peak amplitude propagation 
velocity of exhaust gas pressure waves generated after opening 
successive outlet valves of engine cylinders. During next phase 
of the diagnostic reasoning, harmonic analysis of amplitude 
spectrum of exhaust gas pressure pulses in turbocompressor 
supply channel, is performed by means of the fast Fourier 
transform method [9,10,11]. 

For instance, from numerical data given in the form of 
time - frequency characteristics it can be observed which 
consequences result from elimination of one cylinder from 
operation: drop of the propagation velocity of exhaust gas 
pressure waves in turbocompressor supply channel as well as 
the unwanted concentration of amplitude spectrum of exhaust 
gas pressure pulses in the channel. In the amplitude spectrum, 
the fundamental frequency amplitude dominates that speaks 
about a significant unbalance of gas forces in engine cylinders, 
which makes destructive action on engine construction of 
unbalanced inertia forces generated by rotating and to-and-fro 
moving masses, increasing. The defined diagnostic measure in 
the form of the ratio of fundamental frequency amplitude and 
34 harmonic amplitude in vibration spectrum, for which trend 
analysis of changes of values is carried out during process of 
operation, makes it possible to assess technical state of elements 
of the PRC system, as well as of turbocharging system of the 
engine not fitted with standard indicator valves [11]. 


DIAGNOSTICS OF FUEL SUPPLY 
SYSTEMS OF ENGINE 


In the frame of the research project: „Influence of adjustment 
of fuel supply system of ship piston combustion engine on 
spectrum of shaftline torsional vibrations ” an original technical 
design was elaborated of a test stand for energy — orientated 
investigations of ship engines, for which an universal recorder 
of ship propulsion system performance and efficiency, based on 
measurements of fast-varying quantities, was built — Fig. 4 [2, 
23]. In present, laboratory tests of the recorder are under way 
(on Sulzer 6AL20/24 engine), also it is tested in real conditions 
(on Sulzer 6AL25/30 main engine of H-5 tug). The tests cover 
simultaneous measuring and recording shaftline torsional 
vibrations and fuel consumption of both the engines. They are 
aimed at making assessment of influence of engine technical 
state changes on shaftline torsional vibration spectrum as well 
as fuel consumption. The experimental simulation consisted in 
disturbing the steady-state operation of engine by excluding 
its one cylinder from work. This way was simulated engine 
operation in the state of partial technical serviceability. 

The recorded runs of ship shaftline torsional vibrations are 
subjected to spectral analysis with the use of the fast Fourier’s 
transform (FFT) [2]. This way the diagnostic measures are 
determined for assessing technical state of engine cylinder 
systems and fuel supply system. An original set of methods for 
indirect assessing runs of gas-dynamic processes in cylinders, 
applicable also in diagnosing ship engines (and not only) in the 
case when direct measuring the pressure inside cylinders (due 
to lack of indicator valves) is not possible, is now in the final 
stage of elaboration. The set of the methods is scientifically 
investigated in the doctor thesis of M.Sc. Burski, St. 

The achieved results encourage to continue further 
experimental tests on influence of regulation of fuel supply 
system and timing gear as well as changes in design structure 
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Fig. 3. The stand for the measuring of gas pressure pulses in exhaust gas system — results of diagnostic tests of turbocharging system 


of turbine nozzle (i.e. simulation of contamination of blade 
passages) on the defined diagnostic measures. 


DIAGNOSTICS OF BEARING SYSTEMS 


In the Institute is also under way an intensive research on 
the methods for technical state assessment of engine bearing 
system, based on the measuring of pressure inside cylinders 
under idle run. The research is aimed at finding a relationship 
between technical state changes of bearing system of piston of 
ship combustion engine (leading to its increased mechanical 
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losses), is based on indication of idle running engine at 
uncoupled driving unit. 

Passive experimental investigations on the problem have 
been carried out in the frame of the doctor thesis of M.Sc. 
Wontk, L., since 1997 [23]. In present, in order to verify the 
proposed set of methods for diagnosing ship engines, active 
experiments on the laboratory engine where real defects of 
main and crankshaft bearings have been introduced, are under 
way. For realization of that phase of the research are used sets 
of bearing sleeves having wear of different kind and degree, 
dismounted from real objects out of service — Fig. 5. 
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Fig. 4. The stand for measuring torsional vibrations of ship shaftlines — results of the indirect testing of runs of processes occurring inside cylinders 


As the realized research results indicate, the idle- run 
indicated power constitutes a sensitive diagnostic parameters 
of engine. However it does not allow to unambigously localize 
sources of increased mechanical losses but only to signalize 
their possible occurence. Hence in order to unambigously 
determine whether the changeable technical state of bearings 
is responsible for the increased mechanical losses it is 
necessary to verify the diagnostic tests by means of spectral 
tests of metallic particle content in lubricating oil. In present, 
a programmme of preliminary tests aimed at proper choice 
of analytical models of copper and lead content is under 
elaboration. 


DIAGNOSTICS OF TURBINE ENGINES — 
A HISTORICAL OUTLINE 


Since 1986 in the then Institute of Ship Construction 
and Operation, Polish Naval Academy, several scientific 
research projects on diagnosing ship gas turbine engines have 
been realized [3, 23]. Their scope of themes was inspired 
by the efforts leading to elaboration of effective diagnostic 
methods applicable in real service conditions on board a ship. 
A strategic aim was to pass to the operation process of engines 
in accordance with their technical state. As a result the “ Basic 
diagnostic system” was implemented into the operational 
system of Polish Navy ships. Its concept and realization resulted 
from the following projects: 


> The task project no. 148-23/S-SO/93: „Diagnostics of ship 
gas turbine engines”, 
> The research projects: 


> no.9T12D 008 11: „Technical state assessment of ship 
gas turbine engines on the basis of research on their 
dynamic features”; 


> no.0T00A 009 16:,, Assessment methods of co-axiality 
of torque transmission elements of ship propulsion 
systems” ; 

> no. 0 TO0A 062 19: „Elaboration of a diagnostic system 
for ship gas turbine engines installed on the frigates 
being incorporated into Polish Navy”, 

= and also a dozen of tens of R&D projects and expertises 
financially supported by the State Committee for 
Scientific Research and Polish Navy. 


The system which is continuously developed and 
modernized, contains several diagnostic stands which make it 
possible to perform a complex assessment of ship gas turbine 
engine (practically of an arbitrary configuration) as well as 
to elaborate a prediction of its further operation process. The 
system is equipped with special controlling and measuring 
instruments as well as analytical tools, which make it possible 
to perform diagnostic tests of engines during stopovers and sea 
trials of ships. [4,12,13]. 

In this chapter are presented short characteristics of 
experimental test techniques applied to assessment of 
technical state of General Electric LM-2500 engines and 
Zoria DR76 and DR77 engines, presently subjected to the 
diagnostic program in question [24]. Some metrological and 
organizational aspects of realization of necessary diagnostic 
tests as well as measurements of control parameters of gast 
turbine engine in steady states and transient ones (start-up, 
acceleration and deceleration) are also highlighted. Usefulness 
of particular diagnostic methods regarding their effectiveness 
in identifying engine operational unserviceability, has been 
assessed. In present, investigations are under way to make it 
possible to carry out engine operation process - in accordance 
with technical state - of all turbine engines being in use of 
Polish Navy. Each year, on average 14-16 expertises of engine 
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Fig. 5. The stand for the diagnosing of ship engine bearing system — results of tests of bearing sleeves. 
Pm — power of mechanical losses of engine, PTPC — friction drag power lost in TPC (piston - piston rings - cylinder) system, 
PE — friction drag power lost in bearing system, PW — ventilation drag loss power, 
PZ — friction drag power lost in mechanisms coupled with cranshaft and timing gear. 


technical state are performed, on the basis of diagnostic tests 
of the main functional systems of the engines together with 
their control and safety devices [23]. 


ENDOSCOPIC DIAGNOSTICS 
OF ENGINE FLOW PART 


For endoscopic examinations of gas turbine engines is used 
an extended optical set which makes visual inspection and 
preparation of photographic documentation of internal elements 
of flow part of engine though inspection openings of abt. 5 mm 
in diameter, possible. For carrying out dimensional analysis 
of internal defects of engine elements, their visualization and 


52 POLISH MARITIME RESEARCH, No 2/2008 


documentation in data base a Camedia C—2500L digital camera 
of OLYMPUS firm, is used. The camera is connnected with 
a borescope or fibrescope by means of special links (adapters). 
Fig. 6. presents the endoscopic diagnostic set used in the basic 
diagnostic system of the Institute. 

The borescopes of stiff lense system of various lengths 
make it possible to do observations within side sectors and 
head sector in a wide range of angles of view. The borescopes 
are very comfortable in use during visual inspection of eges 
of guide vanes and rotor blades. In order to overhaul all rotor 
blades the inspection should be connected with simultaneous 
turnining the rotor by hand. Endoscopic examinations of ship 
engines are performed during preventive surveys once a year 


at least, for current assessment of engine technical state, in 
the case when to extend the between-repair time period is 
necessary, as well as in the case of an increased vibration 
level, occurrence of metal particles in lubricating oil, jump 
deviations from trend curves of changes of exhaust gas 
temperature, excessive smoking, control of flow part washing 
effectiveness etc. 


Fig. 6. The diagnostic set of OLYMPUS and STORZ firms: 
1 — set of borescopes, 2 — fibrescope, 3 — set of lighting sources, 
4 — digital camera, 5 — photoprinter 


On the basis of multi-year endoscopic examinations of 


ship engines, the assessment methods of their technical state 
in service conditions were elaborated. They cover a necessary 


Ist stage of gas generator turbine 
of LM-2500 engine — traces of erosion 
on rotor blade surface 


6th stage of propulsion turbine 
of LM-2500 engine- traces of carbon deposit 
on rotor blade surface 


scope and sequence of inspections of internal spaces, which 
make detecting the defects in particular constructional elements 
of engine’s flow part, possible. 

For each type of the engines used on Polish Navy ships 
detail manuals for performing the diagnostic examinations 
with the use of fibrescopes and borescopes, were prepared. In 
Fig. 7 and 8 example places for insertion of endoscope tips into 
flow part of turbine engines subjected to diagnostic survey, are 
presented. Selected defects of ship engines, identified during 
endoscopic examination in operation are presented in Fig. 9; 
they are systematically published in scientific journals, a.o. in 
[4, 8, 14]. 
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Fig. 7. Endoscopic examination of DR 76 and DR 77 gas turbine engines 
—access to their flow parts SNC — Low Pressure Compressor, SWC — High 
Pressure Compressor, KS — Combustion Chamber, TWC — High Pressure 
Turbine, TNC — Low Pressure Turbine, TN — Propulsion Turbine 


1st stage of LP compressor 
of DR77engine — contamination deposit 
on blade surfaces 


9th stage of compressor 
of LM-2500 engine — dent and crack at the edge 
of attack of rotor blade 


16th stage of compressor 
of LM-2500 engine — dent at the edge 
of attack of rotor blade 


Ist stage of HP turbine 
of DR76 engine — loss 
of material of blade plate 


Fig. 9. Defects in ship engines identified during endoscopic examinations 
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Access to elements of flow part 


Fig. 8. Endoscopic examination of LM—2500 
gas turbine engine — access to its flow part 


PARAMETRIC DIAGNOSTICS OF FLOW 
PART OF ENGINES 


From statistical data gathered by ship engine producers as 
well as reliability analysis based on over-twenty-year lasting 
service experience from operation of engines installed on board 
Polish Navy ships it results that this is the flow part of engines 
which constitutes the most troublesome functional system of 
turbine engine and that the dominating unserviceability state 
is caused by different contaminations of blade passages of 
compressors and turbines [3, 4, 7]. During operation process of 
ship engines, is systematically performed assessment of their 
technical state in accordance with trend analysis of changes 
of energy characteristics of their flow parts [12]; this is the 
working medium entalpy flux distribution along flow passage, 
Hex- f (L), prepared for a given load range, usually of 
1.0 P om» (Fig. 10). 


Fig. 10. Distribution of working medium entalpy flux in characteristic 
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control cuts of flow part of engine [7]. 


The so called degree of split of energy characteristics in 
particular control cuts of engine flow part, x-x, determined on 
the basis of measurements of thermal-gas-dynamic parameters 
of working medium, serves as the measure of the generalized 
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diagnostic parameter calculated by means of the following 
relation: 


AH Xx = Hispo = Hx x(0) (1) 


Another method of assessment of technical state of engine 
flow part consists in analysing its dynamic features. As 
a result of degradation of technical state of blade passages of 
compressors and turbines, inertia features of unsteady thermal 
and flow processes in engine are changed - and in consequence 
- also character of gas-dynamic interaction of its rotor units. 
The performed tests prove that a crucial diagnostic problem 
is the experimental determination of changeable capability of 
engine to cumulating and dispersing energy and substances in 
its mechanical and flow systems [13]. The so far performed 
efforts of the Institute’s research team on the problem in 
question resulted in the monograph devoted to diagnostic use of 
acceleration and deceleration processes of multi-shaft engines, 
as well as in promotion of two doctors specialized in diagnostic 
investigations of starting-up processes of ship engines, namely: 
Dr Wioch, J.and Dr Pojawa, B. [13, 19, 22]. 

In present, in the Institute are under way research projects 
aimed at elaboration of a method for analytical determination 
of energy characteristics of engine flow passage, represented 
in the spatial coordinate frame of time and pivot angle of 
guide vanes of some first stages of compressor, considered to 
be the independent variables of occurring unsteady processes: 
Hy yx =f, Tt, a), 

The so arranged characteristics are used for assessment of 
changes of engine dynamic features for diagnostic purposes. 
In this case, an important difficulty is to measure a necessary 
number of control parameters with no delay and with required 
accuracy, as well as the necessity to model characteristics of the 
compressor of variable flow geometry, e.g. by using controlled 
guide vanes of its first stages — Fig. 11. 
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Fig 11. Change of range of operation observed in compressor 
characteristics, caused by interaction of controlled guide vane; 
ABB’B”-— points of cooperation of the compressor with the network 
in steady states; ACC’B”’— points of cooperation of the compressor with 
the network during engines acceleration; B’’D’DA — points of cooperation 
of the compressor with the network during engine 5 deceleration. 


Measurement and recording process of thermal-gas- 
dynamic parameters is carried out with the use of special 
analyzers designed and built in the Institute [23]. 


The problem of model investigations of energy processes 
occurring in axial compressors of turbine engines of variable 
geometry of flow passages, has been considered in the doctor 
thesis of M.Sc. Wirkowski, P. [23]. Complex metrological 
aspects associated first of all with necessity of mounting many 
gauges for measuring pressure and temperature of flowing 
medium, seriously limit possible applications of the method 
to operation process of engine on board a ship. In the engines 
subjected to the diagnostic survey, accessible places for fitting 
additional measuring gauges of low inertia are very limited 
in number. For this reason further work on the problem will 
be carried out by using laboratory stand for testing GTD350 
aircraft engine. The results have been achieved so far encourage 
to carry out experimental research in the area of influence of 
regulation of the fuel supply and starting-up systems of the 
engine, as well as influence of changes in compressor’s design 
structure (by simulation of contamination of blade passages and 
guide vane angle changes) on the diagnostic measures defined 
for engine dynamic behaviour. 


RESEARCH ON VIBRATIONS 


Turbine engine vibrations result from loss of stability of 
its mechanical system, which means that resonance of global 
character appears, i.e. the whole engine vibrates, - in this case 
it is a designer’s problem, or that of local character when one 
of the fundamental vibration frequencies of engine subsystem, 
f, approches the frequencies of periodically changeable forces 
exciting vibrations, La 

The loss of stability of the mechanical system results from 
the following unfavourable phenomena which can occur and 
develop in the engine: the slow worsening of unbalance state 
of rotor systems caused by sedimentation of contaminations, 
erosion and corrosion of blade passages, as well as bent shafts, 
sudden increasing of unbalance of rotor systems caused by 
loss of a part of a rotating element, e.g. blade, exceedance 
of allowable load exerted on bearings, occurrence of friction 
between rotors and engine casing elements in labyrinth seals, 
occurrence of self-excited vibrations due to pulsatory flow of 
working medium (resulting from non-uniform distribution 
of temperature and flow velocity along flow passage 
circumference). 

The phenomenon of resonance is of a great impact on 
reliability and durability of engine as it results in material 
fatigue and consequently — cracks in its structural elements. 
There are two ways leading to elimination of resonance 
vibrations: 

0 Making constructional improvements in order to change 
fundamental frequency of vibrations or that of exciting forces 
(a change of stiffness and mass of rotating elements), 

Removal - from engine operation range — of that rotational 
speed at which the resonance occurs. 

In ship engine operation the basic method to eliminate 
vibrations is the systematic cleansing (washing) of its flow 
part and adjusting of gas temperature field in the control cut 
behind the combustion chamber (usually — behind gas generator 
turbines). In Fig. 12 the example results of measurements of 
transverse vibrations of DR76 turbine engine are presented 
together with results of harmonic analysis of excitations from 
the side of the engine’s rotors. 

The research on vibration of gas turbine engines has 
constituted the subject of scientific considerations of Dr. 
Grzadziela, A., who has extended range of his interest onto 
vibrations generated by any arbitrary structural node of ship 
propulsion system of any arbitrary configuration including 
also the systems driven by piston combustion engines. The 


passive experimental research on the problem, which has been 
continuously carried out since 1997, is aimed at permanent 
improving the methods of diagnostic use of acoustic vibrational 
signal for assessment of state of co-axiality of engines and 
remaining elements of ship propulsion system and power 
consumer [5]. 
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Fig. 12. Linear run and harmonic analysis of vibration signal [23] 
RESEARCH ON BEARING SYSTEMS 


Lubrication conditions of bearing units of gas turbine engine 
decide on its reliability and durability. From service experience 
of aircraft engines it results that most serious catastrophes of 
aircrafts are caused by primary, not detected in a proper time, 
defects of shafts and their rolling bearings. Especially, multi- 
shaft constructions of a complex kinematic system show the 
greatest unreliability of the bearing system, resulting from an 
insufficient fatigue strength. 

In heavy waves, due to detrimental action of cyclic loads 
(impacts and shocks), the rolling bearings of ship gas turbine 
engines are exposed to failures. From available literature 
sources it results that most of the occurred defects was always 
preceded by accelerated wear of the most loaded elements of 
friction pairs, i.e. raceways of outer and inner rings, rolling 
elements in the form of balls and rollers as well as distance rings 
and sleeves and labyrinth seals. Such situation is demonstrated 
by increased temperature of bearing’s structural elements, 
elevated noisiness during operation, changes in physial and 
chemical features of lubricating oil, as well as occurrence of 
metallic particles in it [16]. 

Hence the lubricating oil parameters can be considered 
carriers of diagnostic information about technical state of oil- 
lubricated friction elements of gas turbine engine. The necessary 
information can be obtained from results of measurements of 
the following items: 

° Temperature and pressure in lubricating system (lubricating 
oil rate of flow) 

° Lubricating oil consumption 

° Physical and chemical features of of lubricating oil 

(viscosity, flash point, acid and base numbers, water 

content) 

e Concentration of metallic contamination in the oil. 


The tests make it possible to control wear process of 
a given friction node and to localize — by means of quantitative 
and qualitative analysis - its friction-dependent elements (such 
as the above mentioned) during gas turbine engine service. The 
following means and methods of diagnostic action are used: 
+ Magnetic stoppers and electronic devices sygnalling 
metallic filings, mounted in lubricating oil outlet piping 
— possible detection of large ferromagnetic particles up to 
1000 u,m in size 
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+ Method of radioisotope fluorescence — chemical 
composition of metallic contaminations (small particles up 
to 30 um in size) is analyzed on the basis of excitation and 
measurement of intensity of a characteristic X-radiation of 
a given sample. Leading producers of gas turbines apply 
specific codes to indicate the sensitive structural elements 
covered by thin layer of noble metals, that greatly increases 
engine susceptibility to control in service 

+ Optical methods —optical and electronic microscopes used 
for qualitative assessment, providing sample enlargement 
of /60-times (and more) 

+ Detection of water content and metallic impurities of 2,5 
—30um in size. 


The mechanical tests on contaminations (apart from metallic 
also graphite ones) contained in lubricating oil of ship engines 
and reduction gears, have been carried out by the diagnostic 
team of the Institute since 1992. 

The example results of the tests on Fe and Cu content in 
lubricating oil in function of the operation time t of LM-2500 
engine are presented in Fig. 13. 
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Fig. 13. Runs of changes in Fe and Cu contamination of lubricating oil 
in function of operation time of LM-2500 engine, WO — oil replacement [13] 


On the basis of multi-year tests carried out both by 
the Institute and the engines’ producer, as well as relevant 
experience of aeronautics, the limits of operational tolerance 
field were established: for Fe concentration — 6.0 ppm, for Cu 
concentration — 9.0 ppm. Also, trend functions which make it 
possible to predict friction wear intensity process in engines and 
reduction gears, were elaborated. With their use the scheduling 
of lubricating oil replacement and maintenance intervals is 
made. As a result of the reseach work has been carried out so 
far by the diagnostic team the operational manual of diagnostic 
tests of bearing systems for all the types of propulsion engines 
and reduction gears used on board Polish Navy ships, was 
elaborated; the achievement is based on the doctor thesis of 
Dr Mironiuk W. [16]. 


RESEARCH ON FUEL OIL INJECTORS 


The basic task of fuel supply system of ship gas turbine 
engine is to ensure — in all possible load ranges — an adequate 
(qualitatively and quantitatively) fuel flux delivered to 
combustion chamber. Irrespective of an applied design solution, 
fuel oil injectors take crucial place in the system regarding 
reliability considerations. They ensure a high effectiveness 
of fuel combustion process after sufficient spraying the fuel 
in advance. As a result, a fuel aerosol consisted of drops of 
10-20 um in diameter, is produced [1]. Such atomization of 
fuel makes its whole vaporization and effective combustion 
during flow through engine combustion chamber, possible. The 
greater dispersion of drops the shorter time interval necessary 
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for vaporization and production of uniform fuel — air mixture, 
that determines inertance of combustion process. 

In engine operation process in shipboard conditions certain 
disturbances can occur in correct work of particular injectors. 
They result from continuous process of contamination and wear 
of flow passages as well as deformations of characteristics of 
ageing control elements [23,24]. As a consequence, flow rate, 
quality of fuel atomization and geometry of its outflow from 
injector also become changed. Unfavourable deformations 
occur of gas flux temperature field behind combustion chamber, 
which can lead to a local overheating of its construction. Also, 
turbine subsystems, especially their guide vanes being under 
large thermal-gas-dynamic loads of pulsatory flow of working 
medium, are exposed to failures. 

Control of technical state of injectors is necessary in 
order to prevent possible engine failures and the lowering 
of operational effectiveness of injectors. The most popular 
become disassembling-free control methods based on the 
measuring of selected parameters of thermal cycle of engine, 
possible to be taken in shipboard conditions. Fuel flow rate 
measurements, apart from the measurements of selected 
thermal-gas-dynamic parameters characterizing engine load, 
can be also a valuable source of diagnostic information on 
technical state of injectors [1]. 

In ship gas turbine engines the two-passage swirl 
injectors are commonly used, Fig. 14. They ensure a high 
effectiveness of fuel atomization over whole range of engine 
load changes. 
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Fig. 14. Schematic diagram of fuel flow through two -passage swirl 
injector [1]: a — fuel distributor, b - spring, c — Ist passage of injector, 
d-— 2nd passage of injector, 1 and 2 — control cuts 


The below given formula (2) for determination of mass flow 
rate of the fuel delivered to gas turbine engine injector is the 
most often used way of mathematical description of its work 
in steady states of engine operation [1]: 


k=l 
where: 


k — number of injectors 

Ap = p,-p, — drop of pressure in injectors 

DEM — fuel outflow coefficient 

Ayp — es of spayer nozzle cross-section in fuel outflow 
plane. 


In the process of assessment (calculation) of rate of fuel 
flow into engine combustion chamber under operation, an 
important difficulty constitutes determination of the outflow 
coefficient Lp for particular injectors. It is a function of flow 
character defined by Reynolds number, as well as of a form 
and dimensions of sprayer nozzle [23]. 

For this reason it is usually determined on the basis of results 
of experimental tests. For the analyzed type of injector the value 
of u „= 9.92+0.02, as determined on a laboratory stand of “cool 
combustion chamber”, is provided by its producer [24]. The 
value should be considered only as a preliminary estimation. 
In real operation conditions of injector installed in engine, its 
elements are subjected to direct action of fuel combustion zone 
and they are heated up to high temperature (500-600K). Their 
flow passages become “deformed” and achieve a new, real value 
of the fuel outflow coefficient. Another service problem, as well 
as manufacturing one, is the necessity of selecting individual 
sets of injectors for each of the engine, to match similar values 
of mass flow rate for 1* and 2™ flow passage. Unrepeatability 
of their production results in a substantial scatter of flow rate 
values even up to 5-7% for one set of injectors [24]. 

The diagnostic tests of ship engine injectors are aimed at 
determination of an averaged value of fuel outflow coefficient 
for one set of injectors on the basis of measurements of gas- 
dynamic parameters of working medium, as well as fuel 
combustion. 

The so elaborated loading characteristics represent 
a relation between flow rate and mean fuel outflow coefficient 
in function of changes of the load determined on the basis of 
rotational speed of gas generator rotor coupled with fuel pump 
— Fig. 15. 

0.70 


17800 18800 19300 19800 


Fig. 15. Changes of mass flow rate and mean fuel outflow coefficient 
in function of rotational speed of rotor of HP engine DR76 [1] 
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Comparative analysis of the so elaborated characteristics 
related to the reference values recorded during the initial 
phase of operation makes it possible to formulate a diagnosis 
on technical state of flow passages of set of injectors of an 
engine in service. 


SUMMARY 


O The Institute of Ship Construction and Operation, Polish 
Naval Academy, carries out diagnostic tests of piston and 
turbine engines of all kinds, installed on Polish Navy ships. 
On the basis of currently performed investigations it is 
possible to carry out operation of more than 150 engines 
of the total power output close to 500 MW, in accordance 
with their technical state. 


O The test stands, data base as well as diagnostic programs 
are successively amended and extended regarding workshop 
repair and changes of characteristics of engines during 
service. 


O Until now over 1000 expertises of ship propulsion systems 
fitted with piston and turbine combustion engines have been 
prepared. 


O Presently carried out research projects are focused on 
development and modernization of the basic diagnostics 
system to the extent which will make it possible to carry 
out diagnostic control of ship piston and turbine engines 
installed on board American-origin frigates (fitted with 
Detroit Diesel 16V149 TI engines and General Electric 
LM-2500 engines), which have been lately incorporated 
into Polish Navy. 
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Determination of location of Top Dead Centre 
and compression ratio value on the basis 
of ship engine indicator diagram 


Stanistaw Polanowski 
Polish Naval University 


ABSTRACT 


In the polytropic model of compression process, exponent of polytropic compression 
curve was replaced with a power polynomial in which the piston travel was used as its 
argument. It was shown that 3" order polynomial is optimum one. In the model were taken 
into account the following: design parameters of cylinder, influence of cylinder wear and 
gas blow-by on compression process, error in determining the pistons Top Dead Centre 
(TDC) location, measurement error due to indicator diagram 5 truncation. The presented 
solution of the non-linear model is based on its partial linearization, use of the least 


squares method as well as on application of the optimum determination methods known in the theory of 

experiment. The model makes it possible to determine TDC location on indicator diagram dealing with 

combustion, determination of total compresion ratio, pressure value of indicator diagram truncation, as 

well as determination of maximum values of compression pressure on the diagrams in which self-ignition 
occurs before reaching the TDC. 


Keywords: Top Dead Centre, compression ratio, ship engine indicator diagram 


INTRODUCTION 


Indication of ship engines has been used for regulating 
and diagnosing purposes since the beginning of their use in 
shipbuilding. Mechanical indicators used for that purpose 
make it possible only to determine maximum pressure values 
of combustion and compression and to roughly assess run 
of working processes in cylinders on the basis of developed 
indicator diagrams. 

Appearance, in the 1970s, of electronic pressure analyzers 
often called also MIP (mean indicated pressure) calculators 
and application of computer technique generally have not 
changed situation in the area of gaining information from 
indicator diagrams. 

In the presently offered pressure analyzers no important 
progress has been done as compared with processing capabilities 
of NK-5 pressure analyzers of Autronica, a Norwegian firm, 
which have been the first and most commonly applied 
instruments of the kind in shipbuilding. In many cases, 
especially as far as portable pressure analyzers are concerned, 
one can speak about a regress in comparison with their 
prototype, i.e. NK-5 model. Most of them is not even capable 
of determining the mean indicated pressure. In any of them no 
advanced processing of indicator diagrams can be carried out 
e.g. in the form of analyzing indicator diagram within pressure 
range as well as determining heat emission characteristics 
which can serve as a basic source of diagnostic information. 
Such situation results from met difficulties in indicating ship 
engines in the domain of crankshaft rotation angle and from 
lack of mathematical models applicable to pressure analyzers 


for automatic analyzing the indicator diagrams. However 
a main abstacle is the difficulty in determining TDC location 
on indicator diagrams with a sufficient accuracy, i.e. not lower 
than +0.3°OWK (degree of Crankshaft Rotation Angle). 
The task is especially difficult in the case of measurements 
carried out on indicator valves because of deformations of 
pressure runs, introduced by gas passages and indicator valves 
themselves. 

In the case of low-speed engines when self-ignition occurs 
after reaching the TDC, peaks of pure compression can be often 
observed, such as shown in the indicator diagram for RTAS2 
engine (Fig. 1). 
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Fig. 1. Typical indicator diagrams of ship engines: - for RTA52 engine 
of p, = 1,8 MPa, n = 137 rpm, for A120/24 engine 
of p, = 1,8MPa, n = 750 rpm. 
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The pure decompression runs on the right from the self- 
ignition points were determined by means of extrapolation 
technique. In such case one assumes that the piston’s TDC 
occurs in the maximum compression pressure point or the 
point of the first dervivative’s transition through zero. Location 
of the TDC by using this method may be burdened with large 
values of errors. Determination of the TDC is connected 
with specific difficulties when self-ignition occurs just in the 
TDC or earlier as in the case of medium-speed engines (see 
Fig. 1 — A20/24 engine). Large difficulties in measuring and 
analyzing are introduced by the problem of assessment of 
cylinder tightness, which is usually performed on the basis of 
maximum compression pressure values. If self-ignition events 
occur before reaching the TDC then it is necessary to stop 
delivering the fuel oil to the tested cylinders, that can be allowed 
in service conditions only in special cases because of possible 
loss of tightness of fuel system. Large variability of load is 
often observed, e.g. in electric generating set driving engines, 
that can make obtaining comparable measurement results of 
compression pressure in particular cylinders, difficult. 


MODEL OF COMPRESSION CURVE 
APPROXIMATION 


The simplest model of compression process is polytropic 
one which is usually presented in the following form: 


pV™= const = c (1) 


where: 
V — gas volume 
m — exponent of polytropic curve, m = const. 


The maximum compression pressure values applied in 
contemporary ship engines, exceed 5 MPa, and in the case 
of highly supercharged ones — 10 MPa. For this reason the 
assumption of m = const will usually lead to excessive errors 
in determining the thermodynamic parameters [1]. For the 
highly supercharged engines values of the mentioned errors 
can even excced 5%. 

Having replaced — in the polytropic model — the compression 
curve exponent m = const with the exponent m, = var, finding 
the logarithm and differentiating the model (1), one obtains 
a non-uniform differential equation of the compression process 
in piston engine, e.g. that in the domain of the piston travel 


S, (2): 
oo a gay 
Pas, dS, =” aS 


It should be added that the above presented model is that of 
polytropic process. Its equation does not posess any solutions 
in the domain of elementary functions. The known solutions of 
the differential equations of the compression process [2, 3] are 
of a cognitive character only and thus they are not applicable 
to diagnostic problems. The main limitations are: that TDC 
location is a priori assumed known, lack of association of the 
mathematical model with technical state of piston-cylinder 
system, as well as the not taking into account the disturbances 
introduced by gas passages, indicator valves and measuring 
gauges. 

If the exponent of compression curve cannot be asssumed 
constant and form of its function is unknown then it is natural 
to assume that the exponent m, can be approximated by the 
power polynomial of a- order [7, 8, 10,], e.g. that in the domain 
of the piston travel S : 


+ =0 (2) 


a 
m,=My+m,S,+...+m,S, (3) 
where: 
m, m,, .... m, — polynomial coefficients. 
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The earlier attempts to express the exponent m, as well as 
the thermodynamical model directly in the domain of crankshaft 
rotation angle were abandoned due to a high complicity of the 
model and associated computational difficulties [5, 9, 11]. 

Results of analyses of approximation quality of compression 
curves by using the least squares method, obtained for ship 
engines of various types for the values of a within the range of 
0 +4 show that the least values of sums of squares of deviations 
were obtained for a = 3 (see Tab. 1) [7, 8, 10]. 


Tab. 1. Influence of the order of power polynomial, a, on percentage 
increase of the sums of squares of deviations as compared with the 
approximation based on a =3; a ,— approximation interval 


Type 
of 
engine 


20| 70-168 
2940 | 840 
3 | SRTAS2 [13900] 1340 
| 4 | 6RTA58 | 6820 | 790 | 3 [0 | 3 | 80-175 


The above given results justified to assume the 3" order 
polynomial to be the model of m, exponent in the domain of 
piston travel [8, 10]. 

The instantaneous volume of gas within cylinder, V, which 
appears in Eq. (1), can be expressed as follows: 


V=V -V tV. tV, tV (4) 
where: 
V, — piston displacement volume, 
V_ — piston volume corrresponding with the piston travel 


done starting from the BDC (Bottom Dead Centre), 

V_ — geometrical volume of compression chamber, 

V, — change of cylinder volume due to wear of elements and 
influence of assembling process, 

V _ — apparent change of cylinder volume due to blow-by of 
gases. 


The quantity V, constitutes the volume by which the 
cylinder volume V should be increased to obtain changes of gas 
parameters resulting from gas losses due to blow-by effects. It 
is a function of piston travel. 

By dividing Eq. (4) by the displacement volume V, 
the particular volumes in Eq. (4) can be transformed to the 
following dimensionless form [8]: 


v=l-=-v tV, tV, t Vy (5) 


The value v, is equal to the dimensionless piston travel s, 
done starting from the BDC, which is expressed by the known 
formula for central crank system: 


s= 40-A" yA ”- sin'a —cosa) (6) 


where: 
a — craknshaft rotation angle counted from the BDC 
à — crank radius/connecting-rod- length ratio 


To crank systems of articulated connecting-rods appropriate 
kinematic relations should be applied. 
By introducing the notation: 


vy =l=-v,=1<s, (7) 
and 
No = Y + vt Vox (8) 


the basic model of compression process can be expressed as 
follows: 


PCV + Vox) “= const = c (9) 


for which — according to the made assumptions — the following 
is valid: 


2 3 
M,= Mot M)S,+M,S,+M,8, (10) 


For diagnostic purposes the notion of the total compression 
ratio £, determinable from the below given expession, is 
introduced: 


_lt+v, 


c 
Ve 


(11) 


where: 
V= Va 6,71) 


The volume v, characterizes influence of gas blow-by and 
changes of compresion chamber volume on compresion ratio 
value. 

The notion ,,total compression ratio” was attributed to 
the quantity £, to distinguish it from the notion of the real 
compression ratio defined as the ratio of the working space 
volume at the instant of closing the ports or valves during 
compression stroke and the compresion chamber volume. 

In the performed calculations the simple linear relationship 
of apparent increase of cylinder volume due to gas blow-by: 

V,, — 0.01 s , was assumed. Assumption of another form of the 
blow-by function would not make calculations more difficult. 

It should be mentioned that in service conditions indication 
process is realized in the domain of piston travel or time. Due 
to influence of gas passages and indicator valves, locations of 
the dynamic TDC measured on engine and that thermodynamic 
determined on indicator diagram can significantly differ to each 
other. For Sulzer A medium-speed engines the differences reach 
3,5 °OWK at the rotational speed of 750 rpm [6]. 

For this reason the angle axis of indicator diagram, a, 
should be shifted relative to the angle axis of kinematic system, 
a, by the TDC location error [7, 8, 10], that is expressed as 
follows: 


a=, + A0 (12) 


where: 
Aa, — the TDC location error. 


Value of the TDC location error is searched for. 

The gauges used for indication of ship engines are generally 
dynamic ones in which constant component is eliminated. The 
effect of truncation of indicator diagrams from below occurs, 
which can be clearly observed on the runs of initial compression 
process in cylinders of SRTAS2 engine (Fig. 2). In this case 
measurements were performed by means of a pressure gauge 
of Kistler firm. 

Hence it should be assumed that the real pressure p in the 
compression curve model (9) is the following sum [7, 8, 10]: 


p=p, Ap, (13) 


where: 
p, — recorded pressure run (indicator diagram) 
Ap, — error of indicator diagram truncation (offset). 


APPROXIMATION OF COMPRESSION 
PRESSURE RUN BY MEANS OF GLOBAL 
MODEL AND DETERMINATION OF 
VALUES OF PARAMETERS 


The global model of approximation of compression curve 
constitutes the set of Eqs. (6), (9), (10), (11), (12) i (13). 


The input quantities and parameters of the model are the 
following: 
+ the recorded indicated pressure p, 
+ the crankshaft rotation angle a, (6) 
+ the design (construction) parameter A, (6). 


The searched-for quantities are the following: 
A values of the constants of compression curve exponent 
m,*m, (10), and of the constant c, (9) 
A value of the piston’ s TDC location error Aa,,,(12) 
the total compression ratio £ (11) 
A error (pressure) of indicator curve truncation, Ap,, (13). 


>» 


To determine the above mentioned parameters was used the 


least squares approximation method whose essence consists 
in finding the minimum value of the functional: 
r A 2 
Sp = MIN È (P,- Ban) (14) 
ny 


where: 
Sp -— sum of squares of deviations for approximated run 
MIN- operator of minimum 


p, — corrected indicated pressure 

Pmn ~ Pressure obtained from approximation 

njo = sample number at the beginning of the approximation 
interval a, 

n —sample number at the end of the approximation 


interval a, 


To determine values of the parameters a linearization of 
the model (9) was applied, realized by finding the logarithms, 
twofold application ofthe least squares criterion (14) as well as 
by using the method of experiment theory [4], in order to obtain 
optimum values of the parameters Aa, ¢, and Ap, regarding 
minimization of the sum of squares expressed by the second 
criterion [7, 8, 10]. 

A number of calculations and analyses was performed for 
a few types of medium-speed and slow-speed ship engines, 
from which mutually confirming results were achieved [7, 8, 
10]. The below given presentation of the obtained results is 
based on the analysis of indicator diagrams of SRTA52 engine 
(new one) recorded during sea trials (Fig. 2). 
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Fig. 2. Runs of compression pressure at the beginning of compression 
and close to the TDC, for particular cylinders of S5RTA52 engine 
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It can be observed (Tab. 2) that the determined values of 
Aa, £, and Ap, for particular cylinders are very close to each 
other for both options of determination, and also in the case of 
including two different values of the parameter a, 


Tab. 2. Results of determination of values of the parameters 
Aa,, €, and Ap, for cylinders of SRTAS2 engine (Fig. 2), for a, = 100 °OWK: 
option no. I — simultaneous determination of Aa, € and Ap, 
option no.2 — simultaneous determination of Aa, €, 
Jor the assumed value of Ap, = 0,27 MPa 


Aa OWK 
Cyl. 2 


Cyl. 3 


The observed differences between values of Ap, for 
particular cylinders may awake some doubts (Tab. 2). Results 
of simulation analysis demostrated that in the case in question 
the estimation errors of Ap, values equal to + 20% lead to 
the determination errors of ¢, values not greater than + 2%. 
Therefore it is justified to assume - for calculations — the same 
A p, values for all cylinders, depending on supercharging 
pressure value. 

Analogous results were obtained also from analyses of 
indicator diagrams for engines of other types [8, 10]. 


DETERMINATION OF THE TDC AND 
ASSESSMENT OF APPROXIMATION 
MODEL ADEQUACY 


For analyzing the problem of determination of piston’s TDC 
were used the indicator diagrams where compression pressure 
maximum was observed (Fig. 2, Fig. 3). 

In the case of 6AL20/24 medium-speeed engine the 
diagrams represented runs of pure compression in tested 
cylinder (Fig. 3) [8] because self-ignition occurred in the engine 
before reaching the TDC, under every load. 

For measurement performing in cylinder the value of Aa, 
is small and equal to about 0.1 °OWK, but for measuring on 
indicator valve the value of Aa, is equal to over -0.9 °OWK 
(Fig. 3). It is caused by disturbances from the side of the gas 
passage and indicator valve itself. It should be mentioned 
that on the above presented diagrams (Fig. 3) the angular lag 
between the run on the indicator valve and that in the cylinder 
could not be observed, though it was equal to 2.7°OWK. 

Identical symptoms as those shown in Fig. 3 take place in 
the case of slow-speed engines (Fig. 4). 

For successive cylinders of SRTA52 engine the following 
values of the error Aa, were achieved: -0.35; -0.25; -0.25; 
-0.25, -0.5 °OWK. Analyzing indicator diagrams of engines 
of different types, in some cases one observed Aa, values 
reaching even 1 °OWK. Hence when the TDC location is 
assumed at the compression pressure maximum point or in the 
point of transition through zero of its 1* -order derivative, one 
should take into account possibility of making large errors in 
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determining the mean indicated pressure and characteristics 
of heat emission. 


Cylinder 
p [MPa OWK] 


p [MPa] 


5.27 
177 178 179 180 181 


a [°OWK] 


Indicator valve 
p [MPa OWK] 


p [MPa] 
5.54 


544 Ý 


a [POWK] 


Fig. 3. Comparison of determination results of the TDC location 
on the runs of pure compression in a cylinder and on indicator valve 
of 6AL20/24 medium — speed engine (50% rated load, n = 750 rpm) for 
the approximation interval a, = 100+180 °, and obtained by using the 
approximation model ; p' , p',— the 1” - order derivatives of the runs 
in question. The value of Ap, was assumed with taking into account 
supercharging pressure value. 


p [MPa] p [MPa OWK] 


178 179 180 181 182 183 
a [OWK] 


Fig . 4. Comparison of determination results of the TDC location 
on the indicator diagram of SRTA52 slow-speed engine (100% rated load, 
n = 130.5 rpm, cylinder no.1) for the approximation interval 
a, = 100+180 °OWK. The value of Ap, = 0.27 MPa was assumed 
with taking into account supercharging pressure value. 


The question arises if the TDC locations determined on 
the basis of approximation by means of the total model are 
not erroneous and if the method makes it possible to determine 
the TDC locations in the case when self-ignition occurs before 
reaching the TDC. For adequacy of the elaborated model speak 
the results of the TDC determination for e.g. SRTA52 engine, 


obtained for different values of the ends of the approximation 
inetrval a, > 173 °OWK (Fig. 5). 


Ada, [POWK] 
2 


172 173 174 175 176 177 178 179 180 
a, [°OWK] 
RTA52 
+ Cyll á Cyll 
E Cyll © Cyll 


Fig. 5. Influence of end values of the approximation interval a „on the 
increments of error values, AAa,, for successive cylinders of S5RTA5S2 
engine, at a, = 100 °OWK. 


#cCyll 


The value of AAa,, did not exceed 0,1 °OWK at a „= 173 
°OWK, and at a , = 174 °OWK it did not exceed 0,05 °OWK 
(Fig. 5). It means that prevailing influence on value of the sum 
of squares of deviations has location of the angle axis of the 
kinematic system relative to the angle axis of the indicator 
diagram. Identical results were obtained for engines of other 
types including medium-speed ones [8, 10]. However for the 
values a „< 165 °OWK the values of AAa,, can exceed 1°OWK, 
which is associated with exceedance of the inflection point (i.e. 
zero value of the 2™ order derivative) which — in the case of 
ship engines- is located close to 167°OWK, counting from the 
BDC. Correctness of the method in question is confirmed by 
results of analysis of influence of a „location on determination 
errors of maximum compression pressure values when location 
of the TDC on indicator diagram is known (Fig. 6). 


150 155 160 165 170 175 180 


a [OWK] 
RTA52 AL20/24 

Cyll #Cyll #Cyll E Cyl 

E Cyll ©Cyll ® Valv 


Fig. 6. Influence of value of a „on value of the determination error 
of maximum compression pressure value, Ôp „„» for successive cylinders 
of SRTA52 engines, as well as for pure compression diagrams 
of 6AL20/24 engine, at a , = 100 °OWK. 


As observed in Fig. 6 the errors of prediction of maximum 
compression pressure values, ÒP a» do not exceed 1% even 
for a, = 150 °OWK, and for a, > 160 °OWK they drop to 
0,5 % and take very small values for a, > 160 °OWK, that 
confirms adequacy of the elaborated model and the method of 
determination of parameters. 

Appropriate choice of initial value of the approximation 
interval a, greatly influences obtained results [8, 10]. The 
values of a, should be generally selected in the points where 
valves or ports are already closed and proper compression is 
under way and the diagram is only a little disturbed. For ship 
engines the optimum values of a, are contained within the range 
of 70+100 °OWK depending on type of such engine [8, 10]. 


ASSESSMENT OF DIAGNOSTIC 
USEFULNESS OF THE PARAMETER €, 


The obtained calculation results showed that for 6AL20/24 
engine under various loads, both for measurements performed 
in cylinder or on indicator valve, the determined values of 
g, are close to each other as compared with the maximum 
compression pressure values, p (Tab. 3). 


Pema? 


Tab. 3. Comparison of the values of ¢, and p „„ determined on the basis 
of indicator diagrams of a 6AL 20/24. laboratory engine 


Place of 
measurement 
Cylinder 
Indicator valve 
Cylinder 
Indicator valve 


50 - (pure 


compression) 


100 


In the case of the measurement performed on indicator 
valve the increase of the value of P ma, by about 35 % resulted 
in the increase of the value of £, by 1,4 % only; and in the case 
of the measurement in cylinder — the situation was similar. 
Analyzing the values of s, (Tab. 4), determined on the basis 
of the indicator diagrams of SRTAS2 engine (Fig. 2), one can 
observe that the value of £, for cylinder no. 4 is significantly 
smaller, and for cylinder no. 2 - significantly greater than those 
for the remaining cylinders. 


Tab. 4. Comparison of the deviations of values of the parameters 
e and p „from their mean values, for successive cylinders of SRTAS2 
engine: 6&,— percentage deviations Of €, OP „~ percentage 
deviations of p, 


emax 


Cylinder 
number 
<% . 4.6 4 i 0.4 


ÖP max’ % 


If for assessment of tightness of cylinders of SRTA 52 engine 
one takes into account the values of Pa, (Fig. 2) then in the case 
of cylinder no. 4 the value of òP a =—2.35 % makes it possible 
to consider the cylinder as that having correct parameters of 
compression process in spite of that for the cylinder: ôe =—5.7 
%. In the case of cylinder no. 2 attention should be focused on 
the small value of 8P a =—0.40 %, whereas de, = 4.6 %, having 
opposite sign against that of òP In the case of cylinder no. 
4 a setting error was made and the valve was closed much later 
than in the remaining cylinders, which can be stated from the 
pressure run occurring at the beginning of compression process 
(Fig. 2). It is also possible to demonstrate that the phenomenon 
was not caused by a greater value of the error Ap, observed in 
that case [8, 10]. 

It should be also observed that the error Aa, greatly 
influences the error in determining £, (Tab. 5). 


Tab. 5. Influence of the error Aa, on the error 6é, in the case 
of RTAS2 engine (for cylinder no. 5) 


In view of the so great influence of the error Aa, on 
the determination error of total compression ratio, d¢,, to 
accurately determine the TDC location becomes especially 
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important. In the case of slow-speed engines if measurements 
are performed carefully both in the aspect of pressure and angle 
axis measurements, to determine the TDC locations burdened 
with the error values below +0.05 °OWK, will be possible. 


CONCLUSIONS 


O The presented model and method of approximation of 
indicator diagrams in the compression range make it 
possible to determine the TDC’s location on indicator 
diagrams also in the case if self-ignition occurs before 
reaching the TDC. 


O The so determined TDC is equivalent to that kinematic. 


O The TDC’s location considered identical with the location of 
maximum value of compression pressure or with zero point 
of 1* order derivative of pressure run curve is burdened with 
the error whose value can range from a few decimal parts 
to 1 °OWK. 


O The presented method makes it possible to determine the 
total compression ratio ¢,. This parameter is more suitable 
for assessing cylinder tightness losses than maximum value 
of compression pressure in view of a weak dependence of 
g, values on engine load level. 


O Accuracy of determination of values of the parameter £, 
is greatly influenced by value of the error of the TDC’s 
location. Even the values of the error Aa, equal to 

+0.3 °OWK will produce the errors in determining £, values, 

of the order of +10 %. 


O In the case of slow-speed engines, on the basis of the 
approximation model it can be also possible to determine 
values of the indicator diagram truncation pressure Ap.. 
However it seems that in the case of service measurements 
it would be better to resign from determination of value of 
the parameter in question and to assume - for calculations 
- its value by taking into account supercharging pressure 
value. 
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NOMENCLATURE 
a — crankshaft rotation angle counted from BDC 
a, — approximation interval 
Q — beginning of approximation interval 
«  ~ end of approximation interval 
ös, — error, i.e. percentage deviation of €, 
ÖP max 7 ETOT, i.e. percentage deviation of P aa 
p, — error of indicator diagram truncation (offset) 
Aa, — error in TDC location 
E, — total compression ratio 
p — pressure, given on corrected indicator diagram 
P, — smoothed run of indicator diagram 
Pyemax — Maximum value of compression pressure 
Pa | — run of pressure obtained from approximation by means of 
a model 
p’, .— 1“ order derivative of run curve of p, 


— 1* order derivative of run curve of p, 
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ABSTRACT 


A proposal of probabilistic interpretation of loading process of self-ignition combustion 
engines, was presented. Probabilistic description of loading the engines was proposed with 
taking into account known parameters (indices) of their operation. It was demonstrated that 
instantaneous load of such engines can be considered a random variable. Attention was 
paid to the fact that the load can be taken as a multi-dimensional random variable. Engine 
load changes during its operation were considered as the loading process and presented 
in the form of a multi-dimensional stochastic process whose states are loads considered 


to be random variables. It was demonstrated that the loading process can be taken as the stochastic one 

of asymptotically independent incements, stationary and ergodic. Justification of the above mentioned 

features of the loading process is contained in the presented hypotheses. It was also demonstrated that 

in testing the load process of self-ignition engines stochastic relation between thermal and mechanical 
loading should be taken into account. 


Keywords: hypothesis, load, self-ignition engine, stochastic process 


INTRODUCTION 


Loads belong to the crucial causes of wear both surface 
(linear) and volumetric one, that results in failures of self- 
ignition engines. Especially unfavourable influence on wear of 
the engines is due to excessive thermal loads in the elements 
which form their working spaces (combustion chambers). 
Though mechanical loads of the engines have not so important 
influence on their wear but if they are excessive (especially 
those exerted on main and crankshaft bearings) they can also 
lead to serious failures called breakdowns [3, 7, 8, 15]. 

It is possible to predict such failures in the case of 
application of suitable diagnostic systems so designed as to 
generate full diagnoses, i.e. these capable of predicting their 
technical state. To elaborate such prediction it is necessary 
to determine a.o. loads which can appear during operation 
of engines. In present, loads of combustion engines have 
been generally analyzed deterministically [7, 9, 16, 17, 18]. 
Such analyses based on a probabilistic approach have been 
performed as a rule in a simplified way. However loads of the 
combustion engines should be modeled in time and space as 
random variables attributed to random events in which loads 
of an appropriate value are measured. The fact that a given 
load measurement is a random event, means that to which an 
appropriate probability should be attributed. Moreover the 
loads should be analyzed in successive instants of operation 
of combustion engines. The so considered loads form the 
process of data. Hence it is necessary to elaborate a concept 
of analyzing and assessing such loads in the probabilistic 
aspect with taking into account the fact that the load changes 
successively occurring one by one in the above mentioned 
instants form a stochastic process. However the instants are 
not random variables but parameters of the process. 


In order to elaborate such concept, probabilistic features 
of the loading process of combustion engines should be first 
determined beginning from analysis of states of the process, 
i.e. the engine loads considered to be random variables. 


ENGINE LOADING CONSIDERED 
AS A RANDOM VARIABLE 


In Introduction it was demonstrated that engine load 
values can not be precisely predicted. Therefore the following 
hypothesis H, can be formulated: ,,Engine load is a random 
variable therefore during successively performed measurements 
its values can be predicted only with a certain probability”. 

Qualitative interpretation of combustion engine load in an 
arbitrary instant t can be presented in the following form: 


QO = QW) + Qu (1) 


where: 
QD) = AUO + Q(t) 
and 
Q,(0 — thermal energy (EC) delivered to engine in the 


instant t 

Q.(t) — thermal energy transferred by engine elements during 
its operation (thermal loading) in the instant t 

Q,,(0 — mechanical energy associated with presence of gravity 
and formation of gas forces, friction and inertia forces 
(mechanical loading) in the instant t 

AU(t)— internal material energy increments, in the instant t, of 
elements through which a part of the thermal energy 
Q(t), denoted Q,(t), penetrates 

Q,(0 — energy emitted to environment in the instant t, by walls 
of engine elements, e.g. those forming working spaces 
(combustion chambers) and other. 
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The above mentioned engine loads can be considered 
random variables.They result first from transformation of 
chemical energy contained in fuel into thermal energy, and next 
- the latter into mechanical one (Fig.1). In the interpretation of 
the transformations it was taken into account that the heat is 
a form of transformation of chemical energy into thermal one, 
and the work — that of thermal energy into mechanical one. 


Ss: Sr 


Fig. 1. Example schematic diagram of energy transformation in self-ignition 
engine: ECh — chemical energy; EC — thermal energy; EM — mechanical 
energy; ET — thermal and mechanical energy losses; SoZS — self-ignition 
engine; OEM — mechanical energy consumer (e.g. ship propeller, electric 

generator, compressor, pump); Q (t) — thermal load in the instant t; 
O,(0 — mechanical load in the instant t. 


In the case when only one parameter (one random quantity) 
is used for load description then the load (value of the loading 
process) is one-dimensional random variable. And, if many 
parameters are used to describe the load (i.e. value of loading 
process) then it is a multi-dimensional variable. The load 
is then described by a set of many (generally n in number) 
random variables. In such case it can be considered to be 
a n-dimensional random function, i.e. a set of many (n) real 
functions which attribute unambigously many numerical values 
to each random event, which load measurement constitutes too 
[1, 3, 4, 10, 14]. 

In service practice the self-ignition engine load is usually 
defined by means of the following parameters: p „~ Maximum 
combustion pressure, t „~ maximum combustion temperature, 
p, — mean effective pressure, c, — mean piston speed, AQ... 
— mean pressure rate velocity. The specified parameters can 
be considered to be values of random variables which can be 
denoted respectively as follows: Pao Thao Po Cp and A® 

Ifthe load is determined by means of such random variables 
as: maximum pressure (P „„) and maximum temperature 
(T aax) the two random variables (P ao Tna) can be considered 
simultaneously. The variables can be taken stepwise. Then the 
quantities Pmax and t, max are arbitrary realizations of the variables. 
Hence the pair (p,,t, ) 1s realization of the two-dimensional 

imax. t max 
random variable (P a> Tna). Simultaneous occurrence of 
the events: Pax = Pim nd T, = t. is determined by the 
probability IP aot ae) In this case all the values Pima 200 tg 
which can occur, must follow the relation [1, 10]: 


y 3 OD imine re =1 (2) 


i=l j=l 
The set of the probabilities P(P; nao tna) constitiutes a two- 
dimensional distribution of the random variable (P 
The probability p(p,,,,,) of the random event P= De 
without taking into account value of the random vapa Tiao 
is equal to the sum of the probabilities p(p ), which 
contains all possible values te therefore: 


q (Pimax) = » q (Pimax tma (3) 
j=l 
The set of the probabilities q(p,_,.) determined according 
to Eq. (3), is a boundary distribution of the random variable 


ax? ae 3 


imax? ov 


max’ 


In practice it A neccesary to determine the conditional 
probability q(p. ) of the random event P= Pima, under 


el tax 
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the condition (assumption): T nax = tima: This implies that an 
engine is most loaded mechanically and thermally when in its 
working spaces (cylinders) maximum pressure and temperature 
values occur. The probability can be determined by using the 
formula: 


q (Pimax Uima 
q(Pimax/tjmax) ey (4) 
q(t mad 
The set of the conditional probabilities APimas/ ta) under 
the same condition (assumption): T =t._, is a conditional 


max “jmax? 


distribution of the random variable P_ under the condition: 
=t.. The sum of the conditional probabilities qẸp; 


max jmax’ 
containing all the possible values p. _ equals one, i.e.: 


imax Mad 


imax 


y (Pimax/t max) =1 (5) 
i=l 


In practice it may so happen that the random variables 
P aax and T pa are independent, e.g. as a result of an incorrectly 


performed "regulation of engine. Then the following relations 
are valid: 


q(p ma “nad = q(p ina (6) 
A(t max Pima) = Q(t nax) (7) 


By taking into account Eq. (6) and (7) in Eq. (4) the 
following is obtained: 


q(p max timad) = UP imax) (trax) (8) 


The random variables P a Tax which satisfy the condition 
(8), are stochastically independent. In the case if the condition 
is not satisfied the random variables P a T na, are stochastically 
dependent (correlated) random variables.As results from Eq. 
(8), in investigations should be taken into account the fact that 
in the case of the independent random variables P ao T nax their 
conditional distributions do not differ from their boundary 
distributions. 

The load described by three parameters, for instance the 
maximum pressure P „» maximum temperature T and mean 
effective pressure p, can be considered in the same way. Then 
the three stepwise random variables Cm T may? Po) Should be 
considered. The three quantities: P ame tmas Md Perma Constitute 
realizations of the random variables. Hence the three (Pima t oe 
p.) is realization of the three-dimensional random variable. 
The probability of simultaneous occurrence of the events: 


P =t ax and P,=p,,, constitutes the probability 


max =p imax? ~ max jmax 


POR rae t ‘jmax? Pad- 
In the case of the taking into account of all the values: p 


t aax and p,, which can occur, the following is valid: 


>>) AOimwtis Pek) =1 (9) 


i=l j=l k=l 

The set of the probabilities q(p,_,., t ‘ong Pa) is 

a three-dimensional distribution of the random variable 
(Prao T max’ P.). 

The probabilities q(p,,,.) of the random variable: 

max = Pima Without taking into account value of the random 

variable T as well as the random variable P., is equal to the 

sum of the” probabilities (Pin t mai? Pa) which covers all the 


possible values of ae and p,,, hence: 


imax? 


Dives) = eS A (Pimge mas Pek) (10) 


jel k=l 


The set of the probabilities q(p,_,.) which follow Eq. (10), 
is a boundary distribution of the random variable P, 


max” 


The presented proposal of description of the load of 
combustion engines is important as the load impacts wear of their 
elements. For this reason, the wear in a given instant of operation 
time of the engines is also a random variable, and when analyzed 
in successive instants of the operation, their values should be 
considered as those following a random process. 


ENGINE LOADING CONSIDERED 
AS A STOCHASTIC PROCESS 


As commonly known, in trying to maintain engine’s 
continuous power rate or its rated power, in successive instants 
engine’s mechanical load values are rather close to each other, 
depending on operational conditions. And, engine’s thermal 
load values can, in presence of the mechanical loads, vary and 
even exceed to a great extent their limit values. It results from 
the dependence of thermal load not only on mechanical one 
but also on a state of engine’s cooling system, lubricating oil 
quality, pressure in lubricating system of bearings, oil quantity 
delivered to lubricate cylinder liners, etc [4, 7, 8, 12]. 

The loading (mechanical and thermal) of engines can be 
considered in a dynamic (short) time interval measured in 
[ms] (i.e. that in which one or at most a few thermodynamic 
cycles of engine occur), or in that quasi-static (long) which is 
equivalent to a correct engine operation time interval usually 
measured in [h]. 

By comparing the loads which occur during successive 
thermodynamic cycles, i.e. during the short time interval 
(t,), one can state that they are different because different 
causes [7, 8, 16, 17, 18]. Along with increasing values of the 
engine operation time t, the differences are also increasing. 
Hence at a given value of the dynamical time interval and in 
investigating the loads during particular cycles in the instants 
ti» to- One can obtain different engine load values. The 
event of obtaining a given value of the load is random one. 
It means that to each instant of the time t, a random variable 
equivalent to the engine load in that instant, can be attributed. 
In the same way to each instant of the quasi-static time (t,) 
a random variable equivalent to the engine load in that instant, 
can be attributed. 

The engine load taken as a random variable in a given 
instant t, constitutes value of the loading process. In turn, the 
process constitutes a run of consecutive load changes causually 
inter-connected in function of time, and possible to be described 
by different random parameters (quantities) whose values can 
be predicted only with a certain probability. Hence the loading 
process is stochastic, i.e. a function whose values for a given 
time instant of engine operation, t, are randomly variable. 
Therefore the following hypothesis H, can be formulated: 
„The engine loading process is stochastic because values of the 
engine loading in a given time instant (which is not a random 
variable) are randomly variable”. 

As results from the theory of stochastic processes the set of 
such instants is that of parameters of the process [6, 14]. 

The engine loading process is stochastic one of continuous 
realizations. However by its discretization the stochastic 
process of constant realizations within particular time intervals 
of engine operation, can be achieved. The so modified loading 
process (process of load changes) can be described in the form 
of a semi- Markovian model which has a strictly determined 
set of loading states, [4, 6, 14]. 

The engine loading can be characterized by means of 
various quantities (parameters, indices). Generally, it can be 
expressed in the following form (11): 


QO = QUO, QO] 


engine’s loading 


(11) 


where: 


Q = 


engine’s mechanical loading 
engine’s thermal loading 
— engine’s operation time interval. 


Qu 
Q: 
t 


In empirical loading investigations at least two stochastic 
processes: {Q (t): t > O}and {Q_(t): t > O}can be considered. 
They are components of the vectorial process {Q(t: t > 0}, 
[4, 9, 16, 17, 18]. 

In the case of the loads Q„ and Q, it can be possible to 
write that they are vectors of the following components in 
a given instant t: 


Qu : Pies P2 Pes Ces O: Op Ps Poo] (12) 
Qelqa V TOPe es Tra Ta Ty Qed (13) 
where: 
Paa ~ Maximum combustion pressure; p, — combustion 


pressure; p, — mean effective pressure (p, = NP; 1, — 
mechanical efficiency; p, — mean indicated pressure) 
c, — mean piston speed; @ - degree of isochoric pressure 
increment; P, instantaneous pressure increase rate; 
n-— engine crankshaft rotational speed; P_— gas pressure 
force; P, — inertia force; q — thermal flux (energy) 
density; VT — temperature gradient; p — degree of initial 
isobaric decompression; T „a-~ maximum combustion 
temperature; T, — combustion pressure; T , — exhaust 
gas temperature; Q — heat flux; T, — oil temperature; 
T — cooling water temperature. 


As results from Eqs. (11) + (13), the engine loading depends 
on many quantities (parameters, indices), hence it can be 
considered a stochastic process which contains superposition 
(composition) of particular individual processes (in the simplest 
case — that of Q„ and Q,). Such interpretation of the loading 
is suitable for general considerations but it may be rather not 
useful for practical purposes. 

The loading can be also considered to be a vectorial process 
whose components are random variables which characterize 
load values in particular instants of engine operation time, 
described in the form of Eqs. (12) and (13). 

Generally, the loading can be understood as a process 
whose states can be taken into account in the form of random 
variables. The process, as observed from Eqs. (12) and (13), 
is mult-dimesional because it can be determined only when 
load parameters which form a complete set of parameters, 
are measured and probabilities of obtaining their values - 
determined. 

The complete set of load parameters (indices) can be 
interpreted as that which contains all the parameters necessary 
to determine engine’s loading. 

In every instant t,,(1 = 1, 2, ....n), of the dynamic time 
interval t,, n following random variables: Q, ,Q, ,...,Q, , can 
be considered. The variables have been fui'ther“denoted: O; 
Q>- --» Q, They can stand for n different quantities (parameters, 
indices) of investigated loading (e.g. Q, = Prao Q= Tins Q3 
=p, Q, =q, etc). 

If particular load features are simultaneously considered (the 
load can be characterized by an excessive combustion pressure, 
excessive combustion temperature, excessive heat flux, etc) 
the set of n random variables, being a n - dimensional random 
variable, is obtained. Hence at a distinguished instant t the 
loading Q(t) can be considered to be a n- dimensional random 
variable Q, (further marked Q). Therefore the so considered 
loading can be determined by a set of n real unambiguous 
functions which attribute numerical values to every random 
event (i.e. occurence of a given load value). 
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Generally, in order to simplify further considerations, 
they can be limited by assuming that the random variable Q 
will be a two-dimensional random function (Q,,, Qo). For the 
reason that measurements are periodically executed it can be 
assumed that the random variables Q„ and Q, are stepwise 
(non-continuous) random ones. Realizations of the variables 
are respectively the quantities: q,,, and q,,. Hence the pair 
(du; q,,) is realization of the random variable (Q-Q): The 
variable takes the values (q,,,, q,,) with a determined probability 
Plame Ic) Which is that of simultaneous occurrence of the 
events: 


Qu = Ini ANd Qe = de, 


The set of the above mentioned probabilities p(q,,. qc) 
constitutes the two-dimensional distribution of the random 
variable (Q,,, Qo). In a similar way as in the preceding 
considerations which have made it possible to formulate Eqs. 
(2) and (3), the mentioned distribution Plau; qdo) satisfies the 
following condition: 


22 PCm,-dc,) =! 


i=l j=l 
The boundary distribution of the random variable Q,, is 
as follows: 


(14) 


Plam) => Pm -4c;) (15) 


=l 
and, the distribution of The random variable Q, has the 
following form: 
P(dc;) =>) P(m,-Ic) (16) 
jel 

From the investigations have been performed so far it results 
that some quantities which characterize the loading, e.g. p, €, 
[3, 16, 17, 18] characterize both mechanical and thermal one. 
Therefore it is obvious that there are relations between the 
mechanical and thermal load. For the reason that they constitute 
random processes a stochastic relationship between them should 
be expected. Hence to explain the relationship the following 
hypothesis H, can be formulated: „Between the mechanical 
load Q,,(t) and the thermal load Q (t) a stochastic relationship 
occurs because determined variants of one of the variables are 
accompanied by different variants of the other.” 

Therefore it yields that the relationship between the loads: 
(Q,,(0) and Q_(t)) cannot be described by using the common 
method of algebraic equations. This seems true as the loading 
depends on a large number of factors including those non- 
measurable [2, 3, 12, 14]: 

During the main energy transformation process in 
combustion engine (Fig. 1) thermal energy is transformed into 
mechanical one, but not inversely. Hence the thermal load (Q,) 
can be conventionally assumed an independent variable, and the 
mechanical load (Q,,) — a dependent variable. Analogously, it 
can be assumed that Qep Qez- -> Qop are independent variables, 
but Quai Que a6 Qu, - dependent ones. 

A degree in which the load Q,, is determined - either by 
the load Q., or Q, G=1,2,...,n) considered to be independent 
variables, may be'very different. In practice it may happen 
that one independent variable (Q,.) almost fully determines the 
dependent variable (Q,,). However it may also happen that a few 
independent variables (Q, ) only to a small extent influence the 
dependent variable (Q,,).’From the above said it results that 
there is a necessity of taking into account an intensity (force) 
of the stochastic relationship between Q, and Q,,. 

The intesity (force) of the stochastic relationship between 
Q,,(t) and Q.(t) can be determined, during empirical 
investigations, by using the relationship [10]: 
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2 
T2 c= Taye l 
M CM NJ(k-Dd-) (17) 
where: 
k — number of variants of the variable Q, 
1 — number of variants of the variable Q, 


N — boundary number of the variable Q„ or Qe 
x — value calculated from chi- square formula 
To — Czuprow’s convergence coefficent. 


It can be proved [10] that T „c takes values from the interval 
[0, 1]. The coefficient equals zero (Tuc = 0) if no relationship 
between values of the process (Q,, and Q.) occurs, and if it 
is equal to 1 (Tuc = 1) a functional relationship takes place. 
Hence from the hypothesis H, on the relationship between the 
loads Q„ and Q, the following consequence K can be derived: 
le Tey? 0= Tye? 1. 

Assuming that the syntactic implication H = K is true, to 
verify the hypothesis the following methods can be applied 
[11, 13]: 
€ the reductive reasoning which proceeds according to the 


scheme: 
(K,H>K)-H (18) 
€ the „modus tollens” rule which proceeds in line with the 
scheme: 
(~K,H=>K)}~H (19) 


where: the symbols |} — sign that the hypothesis H is true. 


In this case, application of the reductive reasoning can be 

justified by the following: 

ve Mechanical load increasing as a rule always results in 
thermal load increasing as an increase of the torque M, 
requires an increase of the fuel charge AG, i.e. that of 
the chemical energy contained in the charge. The charge 
increase makes it possible to develop a greater thermal 
energy inside the engine combustion chamber (working 
space). A part of the energy is transformed into the 
mechanical energy (EM) in the form of work (Fig.1) and 
the remaining part, the lost energy (ET), is transferred by 
the chamber walls to cooling medium and - together with 
exhaust gas - to the environment. 

ve Despite the mechanical load does not undergo large 
changes, the thermal load can be significantly increased 
due to worsening the cooling conditions, e.g. due to 
sedimentation of mineral deposits on cooling space surface, 
air penetration to the space, etc. 


Despite the reductive reasoning does not guarantee that the 
conclusion that q,, has increased is true, the relation: (q, = q,) 
N qe is true. However it does not mean that, if the implication: 
(q,, = qo) on the interpretation that if the mechanical load q,, 
increases the thermal load q,.” also increases, it can not be 
assumed true that if q, increases then also q,, increases, i.e.: 


[au > Ie) N qel => Any (20) 


And, application of the „modus tollens” rule results from 
that the scheme of reasoning on the hypothesis H, is reliable 
and leads to its falsification. However it is so only in the case 
if empirical investigations are truly capable of undeniable 
contradicting the hypothesis H,, i.e. that the mechanical load 
will increase. 

Itis possible to state generally whether the random variables 
Qu and Q, are mutually dependent, by investigating the 
probabilities determined with the use of the general equations: 
(14) + (16) or those detailed: (6) + (8). 


In the case when - as in the detailed considerations — the 
conditional distributions of the variables do not differ from 
their boundary distributions, the following is valid: 


P(am, *Ic,) = PC, ) Pe, ) (21) 


It means that the random variables Q,, and Q, are 
independent. And, if Eq. (21) is not satisfied they are 
stochastically dependent. In this case verification of the 
hypothesis H can follow the schemes similar to those presented 
in the form of Eqs. (18) and (19). 

Investigation of stochastic relationship between the random 
variables Q„ and Q, is difficult because distributions of the 
investigated loads (including conditional distributions) can 
differ in many features (variability, concentration, asymmetry 
etc). Hence it is much easier to examine statistical relationship 
between the mentioned random variables Q, and Q,- 

It should be observed that if no stochastic relationship 
between the random variables Q,, and Q, takes place then any 
statistical relationship does not take place between them too. 
However reverse implication is not true. It directly results from 
that the same conditional mean values are related to the same 
conditional distributions, whereas to the same conditional mean 
values the same distributions are not to be related [1, 10]. 

The statistical relationship is distinguished by that to 
particular mechanical load values to attribute a mean value of 
thermal load, or reversely, is possible. 

Investigation of statistical relationship between any random 
variables is realized by means of statistical ordered series 
or appropriate regression equations which make it possible 
to determine whether a correlation takes place between the 
mentioned loads and - if it does -whether it is positive or 
negative. In practice particular realizations of self-ignition 
engine loads are observed. In the case if realizations of load 
for a given time instant t differ only a little, i.e. if its particular 
values for a given time instant t constitute (for increasing n, 1.e. 
number of load measurements) the series of random variables 
Q, (n= 1,2,...), then the series is stochastically convergent to 
zero under the following condition, [1, 10]: 

limP{|Q, > €} =0 (22) 
for an arbitrary number e > 0. 

The above is equivalent to the statistical convergence of 
the series of random variables {Q_} toward the constant value 
q # 0, hence — to the statistical convergence of the random 
variables {Q, — q} toward zero, that can be expressed as 
follows [1, 10]: 


lim P{|Q, -—q|>e} =0 


lim )Q, = q 


noo 

The relation (24) makes it possible to consider loads as 
a deterministic process, i.e. that statistical having prediction 
error equal to zero. In such process the past fully determines 
the future.The simplification consisting in the assumption that 
“the loading constitutes a deterministic process” is commonly 
applied in investigating the engine loads. However it makes 
load assessment to be fully adequate to reality both as regards 
to its features, especially to be capable of rational forming 
its run, i.e. its rational steering, and - if possible - its optimal 
control, difficult. 

Comparing the loads occurring in successive cycles 
(number of cycles n > 4), hence within the time interval t, 
(short), one can state that they are different. The differences 
are increasing along with increasing n, i.e. when considering 
the load within the time interval t (long). For this reason 


(23) 


or 
(24) 


different load values are obtained at a determined value of 
the time interval t, and during investigation of the load in 
particular cycles at the instants t,,,t,,,...The fact of obtaining 
a given (expected) load value is a random event. This is such 
event because as a result of determination of the same empirical 
conditions the expected load value can occur, but it can also 
not occur. It means that a random variable can be attributed to 
every instant t,. Similarly, to every instant ta relevant random 
variable can be attributed. 

From the above presented considerations it results that the 
following hypothesis H, can be formulated: „The loading is 
a process of asymptotically independent increments because 
along with increasing partition between the time intervals 
within which the load is investigated (load measurements are 
made) its values are less and less mutually dependent’. 

Another feature which characterizes load changes consists 
in that observed load values do not show any monotonic 
changes both in the dynamic time (t,), i.e. that necessary to 
realize one engine cycle and in the quasi-static time (t,) in 
which the engine delivers power to a consumer, e.g. screw 
propeller. Hence it can be assumed that peak values of the 
quantities which characterize the loading, appear at random. 
Therefore to unambigously (precisely) predict an instant of its 
occurrence, is not possible. The lack of monotony feature of 
engine load changeability makes it possible to formulate the 
hypothesis H, as follows: „The loading is a stationary process 
because in a longer time the monotony feature of engine load 
changes is lacking”. 

The load stationarity (in a broader sense) means that in 
each case all multi-dimensional probability density functions 
depend only on the mutual distance of the instants T,, T,,..., T, 
and do not depend on themselves [2, 5, 12, 14]. Hence one- 
dimensional probability density function of load values does 
not depend on the instant which the value corresponds with, 
and two-dimensional one depends only on difference of the 
instants in which observed load values have appeared. And, 
in a narrower sense, the stationary loading (fully stationary) 
is understood as that whose all statistical moments of higher 
orders as well as total moments of the loading (considered as 
a process) does not depend on time. In the case of the stationary 
process (in a narrower sense) the expected value m(t), variance 
V(t), auto-correlation A(t,,t,) and auto-covariance K(t,,t,) do 
not change. And, the stationary process in a broader sense is 
characterized by that m(t) = m = const as well as A(t,,t,) = 
A*(t,-t,) = A*(r). In practice, the stationarity of the loading 
in a broader sense is of importance. 

To confirm the presented features of the stochastic process 
which the engine loading constitutes, it is necessary to perform 
relevant empirical investigations. From the tests of self- 
ignition engines, have been carried out so far, it results that 
their loading is changing continuously so that its particular 
values measured after passing very short time intervals, are 
strongly correlated to each other. However, if the time partition 
between measurements increases, correlation between the 
loads decreases. Hence the loading values measured within 
time intervals (or instants) very distant from each other, can 
be assumed independent. The feature is called the asymptotic 
independence of the load value measured e.g. in the instant 
t,,, and the value measured in the instant t, i.e. when the 
partition (time distance) At = t, ~T, is sufficiently long. The 
so considered asymptotic independence between load values 
measured or calculated in the instants t, and q, , illustrates the 
fact that along with increasing At dependence between them is 
decreasing. And, from the self-ignition engine work principle 
results also that within a longer time of its correct operation 
there is no (and can not be) monotonously increasing or 
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decreasing changes. Hence it can be assumed that maximum 
load values appear at random in determined instants, therefore 
they can be predicted only with a certain probability. This lack 
of monotony feature is called the load stationarity. 

Therefore it can be assumed that verification of the 
presented hypothesis can be performed by checking the 
following consequences which result from it: 


- K; m(t) =m, = const 
K`: A(t.) = A*(t7t,) = A* (At) 


To verify the hypothesis in question the reductive reasoning 
in accordance with Eq. (18) can be used. 

Hence the hypothesis H, of the following content can be 
formed: „The loading of every engine is an ergodic process 
because the loading values do not depend on its initial 
state”. 

The hypothesis is of practical importance as it implies 
that the loading can be cosidered to be a stochastic process in 
which observation probability of the value q(t) belonging to 
the strictly determined set A C R, can be estimated by means 
of the mean time of presence of each of its realizations in 
the set at a long time of observation (measurements). This 
hypothesis - together with the preceding one - explains why 
vibro-acoustic diagnostics of machines in which vibro-acoustic 
signal is assumed stationary and ergodic, is useful in practice. 
As results from the just presented hypotheses, this is an obvious 
assumption as the loading (cause) generates vibro-acoustic 
signal (consequence). 

The statement whether the loading (taken as a stationary 
process) is ergodic, consists in comparing its statistical 
characteristics (expected value, auto-correlation, variance 
etc) achieved as a result of averaging the realization set 
{q,(t)}, with the statistical characteristics obtained as a result 
of averaging a single, sufficiently long realization of the 
loading. 

As results from the hypotheses (from the 3" to the 5%, 
inclusive), the loading investigated in the instants very 
distant from each other, can be considered a purely random 
process, i.e. such stochastic process in which all the random 
variables q(t € R,), in the discrete time (t,), are mutually 
independent. 

The presented opinion on engine load features may lead 
to new possibilities in empirical determining wear- to- load 
relationships. Probabilistic approach to engine load is presented 
in some publications, e.g. [2, 4, 12], where it have been assumed 
that load realization is normal process and its values are random 
variables following the Moivre — Gauss distribution. 


FINAL REMARKS AND CONCLUSIONS 


O In the presented analysis and synthesis of events it was 
demonstrated that the loading of every self-ignition engine 
considered in a given instant of its time of operation (work), 
can be taken as a multi-dimensional random variable. The 
loading analyzed in successive instants of operational time 
of such engines can be considered to be realizations of 
the loading process. Hence the loading process of every 
engine should be investigated under assumption that it is 
a multi-dimensional stochastic process. The hypotheses 
have been proposed why it is possible to assume that the 
loading process of any self-ignition engine can be taken 
as a stochastic process of asymptotically independent 
increments, and stationary and ergodic one, as well as that 
there is stochastic independence between its mechanical 
and thermal loads. Intensity of the stochastic relationship 
can be stated by using the Czuprow’s test of convergence 
[10] during investigations. 
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O For verification of the presented hypotheses the method of 
non-deductive (inductive) reasoning called the reductive 
reasoning, as well as the deductive reasoning method called 
the „modus tollens” rule, have been proposed. 

O To recognize features of the processes it is necessary to form 
relevant mathemetical models by using the system modeling 
and by performing suitable empirical investigations. 

O From the preliminary considerations it results that the 
loading processes of combustion engines are the processes 
whose features can be analyzed by applying the theories of 
decision (controlled) semi-Markov processes [4, 5]. 
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A discrete model of the plate heat exchanger 
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ABSTRACT 


Developing numerical control methods requires precise models of technical machines. The 

article presents a discrete model of a plate heat exchanger. The model was worked out based 

on differential equations to secure high accuracy. Correctness of the model was verified 
using the data for the M10 — MFM cooler produced by Alfa-Laval. 


Keywords: simulation, dynamic model, Z transform, plate cooler, ship facilities 


INTRODUCTION 


Heat exchangers are used in all types of power plants for 
heating or cooling the working medium having the liquid or 
gaseous form. The operation of heat exchangers is crucial for 
keeping accurate temperature of the process. 

In marine power plants three-step main engine cooling 
systems have been introduced, which increased the number 
of coolers and control systems. Plate heat exchangers are in 
common use here. Moreover, the plate construction is also used 
in the designs of condensers, evaporators, and fuel heaters. 

Marine power plant simulators require fast algorithms 
simulating the operation of plate heat exchangers to determine 
new control points for numerous installations. At the same time, 
the simulation is expected to model properly static and dynamic 
properties of the heat transfer process [7, 13, 15]. 

The principles applicable when working out the simulation 
algorithm have been formulated in [10]. The there presented 
procedure is used for working out equations for the fresh 
water/sea water heat exchanger model. 


MODEL OF THE PLATE HEAT 
EXCHANGER 


The plate heat exchanger consists of plates arranged parallel 
to each other, between which the first and second working 
medium flows alternately. The heat transfer process is repeated 
between successive plates. Each plate takes the working 
medium from the main collector and generates repeatable heat 
transfer conditions, including constant medium flow velocity 
and heat transfer conditions represented by the overall heat 
transfer coefficient. In order to develop a mathematical model of 


the plate heat exchanger, a two-plate system with co-current or 
counter-current flow is to be analysed, Fig. 1. When modelling 
the heat transfer, the following simplifying assumptions were 
adopted [2, 3, 4, 12, 14]: 


ve physical parameters of the working media do not depend 
on temperature within the range of temperature changes 
taking place in the exchanger, 

ve the averaged profile of the working medium flow velocity 
is characteristic for the plug flow (turbulent flow), 

ve temperature gradients in y and z directions are equal to 
zero, 

ye thermal capacity of the plate membrane is neglected due to 
its small thickness (g = 0.4 mm). 


Fig. 1. Scheme of the plate heat exchanger segment 
for the counter-current flow 


Taking into account the above assumptions we can write the 
energy balance for an element having length Ax [11, 13]. The 
heat transfer process in the plate heat exchanger is described 
by a system of differential equations (1) and (2): 


oT), y OT) _ ~a[T, (x, t)—T,, (x. 2)] «1 
OT Ox 

aT n(x) Wp Tar) =b, (x,t)-T,, (x, D 
ôr Ox 
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a= (3) 
ZCP, 
b= 2k (4) 


The basic variable is the temperature of the liquid T and T „ with respect to plate length x and time t. The system of equations 
(1) and (2) was solved using the Z transform, and its solution was: 


T, (xa) a 2) +w, daa 27, (x,0) =a-T_(x,z) (5) 


p EE? Tno) = b-Ts(x,z) ie 


= z-l 
Tm(x,z)| b +— |+ 
(2) =| w 


The process is nonlinear due to the flow velocity w and heat transfer conditions k. However, their changes are slower than 
temperature changes in a given discretisation step, and are assumed constant in a single iteration step. This way we arrive at 
a system of first-order differential equations (5) and (6) with constant coefficients. 

To solve the system of equations (5) and (6) we assume a constant non-zero value of the temperature along plate length in 
time t = 0: 


T,(x,0)=T „(x)=const T,.(x,0)=T,,(&)=const (7) 


The general solution of the system of equations (5) and (6) for non-zero initial conditions (7) is a composition of 
functions: 


T; (x, z) = Ae™ + Be™ +T ss (x, z) (8) 
Ta (x,z)= [at kale f Jae” (pa w, r Je” + Tim (x,z) (9) 
aT a aT a` 
where: 
1, =a z+p-q(e+By +y |-2x-4 x (10) 
L Tw, W, 
n =-of z+B- 4+8} +7 |- Il e (11) 
L TW, Win 

= z— tbe az 2(z—1)T,, +zT(bT,, +aT,,,) 

bras 1 Z he ea aa (2 
Tw wain — Tw w,ib (2-1) +(z—-1)T(a+b) 
Z——+az 
Fam (x,2)= bz t4 T,= z(z -1T +zT(bT,, +aT„o) (13) 
Tw, wai Tw.w tt, (z-1) +(z-1)T(a +b) 
and 
= Win =W; 
7 2Tw,,W, an 
p= T2¥m —bw, 1 (15) 
Win “W; 
y= 4abT? > (16) 
(wa =w.) 
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The complete solution will be obtained for known initial conditions. The most typical flow arrangement used in plate coolers 
bases on the counter-flow. The assumed initial conditions include a step temperature change at plate cooler inlet, recorded 
simultaneously on the cooled and cooling medium sides, Fig. 1: 


Ts(x =1,z) = Ae“ +Be™ + Ts (1, z)=Tsi(l, z) (17) 


Taa =0,2)=(1+ 25a, Jaf: s B+ Ton (0,2 = Tm (0.2) (18) 
a a 


aT 


After calculating the integration constants A and B from equations (17) and (18), and placing them to the system of equations 
(8) and (9) we get the solution: 


2 Zl, a. 
T+ ¥ : gee (z+B) en 
= 5 = (z+B+ +B) +7) 
T.2(0,z) =[Ta(,z)-Ta(2)| _ 
1+ id Raids GORIN z- Le 


(z+p+Ve+p) +) 


(19) 
fe Var] vince Sa 
Fabo 2)—TFan( aa) = 2Taw,  z+P+y(z+B) +y z+B+4(z+By + 
—w,) P y cial Aa w 
(z+ B+Ve+py +y) 
iTe (0, z) 
Tm2(1,z) = 
1 a a Vesar] í E U 
+ 
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2 
(z+B+ye+B) +y) 
T(z) 
The original form of functions T., T „ is obtained after applying the definition [6]: 
i an = 1) 
. p 
f, slimo a (21) 
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where: 
at for odd n 
S= 4 
= for even n 
1 -20| 2+8-VleBP + | 
F, (z) = ———— = e 
B+ lz By +7 
ni 3 : k+1Xn-1} bee k 
f T)=0 yee re- gt ain 
aT) ALAI ) ( ) ( ) 2"""kimi(m+k+1\(n—k—1—2my> ý (cy) 
where: 
[> = =! forodd n and even k or evenn and odd k 
S=4 2 
| — : T2 foreven n and even k orodd n and odd k 
F (2) _ l a +] 
(2+B+Ve+B) +7) 
n-2 8 k +2)(n-1) BEF k 
f T) =0.0 “aka eg oom m+1 
se: k x ) ( ) ( ) amimi k+ 2a- k2- 2m]. i (cry) 
2 a forodd n and even k oreven n and odd k 
S= ~ 

= : =2 foreven n and even k orodd n and odd k 

F, (z) 7 y -20| 2+B-y(z+8P +7 | 


(2+B+J@+B) +1) 
a ED peany (oy) 


om 2°"? ktm!(m +k +2}(n-k-2-2m} 


n-2 


f, (nT) = 0,0, 
k= 
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(23) 


(24) 


(25) 


(26) 


After dividing numerically by series we arrive at the response in discrete time instants: 


b 
Tyo (nT) = >) ken-ta oe L Wae W 


—— S 
_ aw,,—bw, 
Wwa =w) tetfar——L yt ft 


Wi Ww, 
Ps É (27) 
f [r-e +f [arh 
W W 
+5 [Ty (2T)= Tyn (0,nT)] m = +T aQ, nT) 
7 | Lt | 
1+f,} nT -——+— os 
Win W, 
2 i | 4 
pare e™ +f | nT +— le™ 
W W 
T,,(0,nT) = $ [T,,(nT)-T,, (nT) H 
In Yy 
” itar- aih ar 
Wa W, (28) 


The excitation functions T, and T „ can be arbitrary, but 
this requires the use of the function convolution operation. For 
the step excitation function, its response is the sum of terms of 
a given series. In case of another function, its approximation 
by a staircase function is recommended. Such approximation 
is used in numerical control algorithms. 

The final solution includes all process variables and the 
time in a discreet form. In one discretisation step the process 
variables take constant values and can be only changed in the 
next discretisation step. Thus we reduce the model to a linear 
form. By modifying parameters in successive steps we arrive as 
a quasi-linear solution to a nonlinear system. If the parameters 
do not have to be changed, a series of consecutive process 
values is calculated. 

Sample results of calculations performed using the above- 
described method are shown in Fig. 2. 


M10-MFM cooler 
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Fig. 2. Interval of the transformed function discretisation 


The simulation of cooling plate operation made use of the 
following input data: 


AW m ~bW | 
f: (nT)- f oT = A + fe WW s 
W W, 


M10 —- MFM cooler, type: fresh water - sea water 


1= 0.675 m 
a= 0.75 [l/s] b = 0.74 [1/5] 
w, =1.23 [m/s] w, = 0.944 [m/s] 


The initial temperature distribution along the plate length: 
T(x) = 39.3°C Tm) = 39.3°C 
The signal assumed at cooler input: 


T,,(t)=75.0°C — T,,(t) =39.3°C 
CONCLUSIONS 


The presented model of the plate heat exchanger reveals 
good agreement with the real results. It can be used in numerical 
simulations, depending on the available simulation time. 

The performed cooler model analyses make the basis for 
formulating the following conclusions: 


O The obtained static and dynamic characteristics are fully 
compatible within the entire range of loads taking place in 
operating practice. 


O The model correctly reacts to changes of working medium 
flow velocities and other parameters. 


O Detailed analysis of the simulated phenomena and processes 
is possible. It can refer to both the control processes, and 
diagnostics of the technical state of the examined cooler. 


O Comparing counter-current and co-current flows produces 


better results, but makes the time in which the steady state 
is reached longer. 
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NOMENCLATURE 


A, B — integration constants 

a,b — heat transfer constant for plate cooler cooled and cooling 
side, respectively [s7] 

c — specific heat of the liquid [J/kgk] 

k  — plate cooler heat transfer coefficient [W/m°K] 

l — plate length [m] 

p  — liquid density [kg/m°] 

T_ — sea water temperature in the plate cooler [K, °C] 

T. — fresh water temperature in the plate cooler [K, °C] 

T  — interval of the transformed function discretisation [s], 


Z(f)=f(2) 53 f(aT)z " 


T — time [s] 

x  — temperature distribution along plate length 
w,, ~ sea water flow velocity [m/s] 
w, — fresh water flow velocity [m/s] 
z  — distance between plates [m] 
Indices 

0 — initial value, 

1 — input boundary value, 

2 — output boundary value, 

m — sea water, cooling medium 

s — fresh water, cooled medium 
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ABSTRACT 


In geological researches several kinds of methods are applied to discovering the natural 
resources. Planes, helicopters and UAVs (UnmAnned Vehicle) are used in researches in large 
areas. The gravity, electromagnetic and magnetic methods, which are used in geological 
researches, are presented in this paper. The special attention was paid to magnetic systems 
installed on mobile platforms. The magnetic field of the Earth obtained from mathematical 
model was compared to the real magnetic field in the selected part of the Baltic Sea. The 
results of the calculations showed that the mathematical model of the Earth ’s magnetic field 


does not consider local magnetic anomalies. The strong local magnetic anomalies cause serious problems 
with detection of underwater objects. Special problems appear in the magnetic system on a helicopter, 
which are presented in this paper. 


Keywords: geological researches, gravity, electromagnetic, magnetic methods 


INTRODUCTION 


In geological researches several kinds of methods are 
applied to discovering the natural resources. Planes, helicopters 
and UAVs (UnmAnned Vehicle) are used in researches in large 
areas. The gravity, electromagnetic and magnetic methods are 
applied to researches. In the gravity method the extremely 
small changes of the Earth’s gravity are measured [1]. The 
sensitivity of contemporary gravimeters is less then 1 mGal 
(1 Gal = 10° m/s’). The gravimeters allow to measure the 
changes of the Earth’s gravity on level close to 10° g. One of 
a few leading producers of extremely sensitive gravimeters 
is Gravimetric Technologies Ltd. a company in Moscow [2]. 
The GT-1A gravimeter installed on the Cessna 404 plane was 
shown in fig.1 [3]. 

Deposits of natural resources are discovered also by 
applying the electromagnetic method. The first electromagnetic 
systems appeared and were developed in Scandinavia, the USA 
and Canada in the 20s of XX century. The electromagnetic 
method is applied to measure the conductivity of the soil. The 
electromagnetic systems are installed on planes or helicopters. 
A big coil is towed by a helicopter or carried by a plane. The 
current impulse in the coil generates strong magnetic field 
(primary field), which penetrates layers of the Earth (fig. 2). 
The time-varying field produces eddy currents in the soil. 
After the switch off of the current in the coil they are only 
the eddy currents that produce the magnetic field (secondary 


Fig. 1. The GT-1A gravimeter installed on Cessna 404 plane [12] 
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field). The time-varying secondary field produces emf voltage 
in the measure coil. The induced emf voltage depends on 
a conductivity of the soil. The analysis of the conductivity of 
the ground allows to detect natural resources. 


ee | 


Primary field 


Fig. 2. The electromagnetic method 


One of the biggest commercial companies (Fugro) has three 
electromagnetic systems for geological researches: GEOTEM 
(fig. 3), MEGATEM and TEMPEST. The biggest MEGATEM 
system is installed on Dash 7 plane and allows to discover 
the natural resources up to 300 m. The value of the magnetic 
moment generated by the system is 2.2:10° Am’. 


Fig. 3. The GEOTEM electromagnetic system on the plane Casa 212 [4] 


The magnetic method is also applied to geological 
researches. Homogeneity of the Earth’s magnetic field is 
disturbed by different properties of the Earth’s layers. This 
method is also applied to detection of ferromagnetic objects 
(UXO — UneXploded Objects, submarines, wrecks, landmines 
and marine mines). Each ferromagnetic object disturbs the 
homogeneity of the Earth’s magnetic field (fig. 4). The measures 
of the magnetic field near the object allow to detect, to localize 
and to identify it [5]. 


Magnetometer = 
3 va N 


Anomaly field 


eooo 


Object's field 


Fig. 4. The disturbance of the homogeneity of the Earth's 
magnetic field near the ferromagnetic object 


The influence of the strong local magnetic anomaly on the 
magnetic measurements, with systems installed on planes and 
helicopters, was presented in this paper. The mathematical 
model of the Earth’s magnetic field was described. The 
magnetic field based on the mathematical model was compared 
to the real magnetic field of the Baltic Sea. The calculations of 
the magnetic field were carried out in MathCad [6]. 
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EARTH MAGNETIC FIELD- 
MATHEMATICAL MODEL 


The magnetic field of the Earth can be described by 
magnetic scalar potential if j = 0: 


H=-V¥ (1) 


The magnetic flux density over the surface of the Earth is 
B = u H. So the equation: 


VB=0 (2) 
comes down to the form: 
VH=0 (3) 


On the basis of the equations (2) and (3) the magnetic scalar 
potential ¥ satisfies the Laplace’s equation: 


VF =0 (4) 
The Laplace’s equation in spherical coordinate system is: 
bid 

a = i o + 

ro or r'sinð ð (5) 
bor 
r’sin’O Ap? 
The solution of the Laplace’s equation (5) is [7]: 
œ =] 

¥(1,0,) = >) >) Amr + Bint) Y;" 0,0) © 


1=0 m=0 


where: A „ and B „~ constants and 


(21+ 1)d-— m)! 


1 
2 : 
Yeo) =| An(l+m)! J Paton n 


where: P, „~ Legendre polynomials [7]. 


Schmidt developed the formulae describing the magnetic 
scalar potential in 1939 [7]. These formulae allow to calculate 
the approximate value of the components of the magnetic flux 
density at any point of the Earth. The magnetic scalar potential 
of the Earth is: 


Y(r,0,ġ)= 
= ay 5 P” (cos Xg cosmo +h sinmọ) 


0 1=1 m=0 


(8) 


where: P™ — partially normalized Schmidt’s functions [7]. 


The Gauss’ coefficients (g," and h™) (8) are obtained by 
using many sets of observations of the components of the 
magnetic field at many points of the Earth. The components 
of the magnetic flux density of the Earth are (for 1 =- 1, 
1 =1,1 =-1): 

y poz ři 


_ low 9 
B= tA (9) 
-| ow 

ve Fs Ob am) 


_ ow 
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(11) 


The modulus of the magnetic flux density in the Cartesian 
coordinate system is (fig. 5): 


2 2 2 
B. =B, +B? +B2, (12) 
The magnetic inclination is: 
I= arctg re (13) 
(ot Bie 
and declination: 
D =arctg (14) 
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Fig. 5. The components of the Earth's magnetic field 


The characteristic components of the Earth’s magnetic 
field were calculated in MathCad for 4-th degree of Gauss 
coefficients (1 = 4). The isodynamic chart of the total magnetic 
flux density was shown in fig. 6. The isoclinic chart of the 
inclination in degrees was shown in fig. 7. 
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Fig. 6. The isodynamic chart of the total magnetic flux density 
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Fig. 7. The isoclinic chart of the inclination in degrees 
The distribution of the Earth’s magnetic field given by the 


mathematical model can be used for example in demagnetisation 
systems of war ships. This distribution has global character and 


does not include local magnetic anomaly. It is the consequence 
of the limited number of observatories at which the components 
of the magnetic field are measured. The distributions of the local 
magnetic anomalies are different and depend on “geological 
history” of the Earth’s layers. The most popular magnetic 
anomalies in the world occur in the surroundings of Kursk. 
The massive layers of rocks containing the iron ore cause these 
anomalies. The values of the anomalies are about 27 nT 350 km 
over the Earth’s surface. The relatively strong local anomalies 
occur above the Baltic Sea. The distribution of the modulus 
of the magnetic flux density obtained from mathematical 
model (8-12) of the Baltic Sea area was shown in fig. 8. The 
distribution of the magnetic field is approximately linear. The 
real distribution of the magnetic field in the same area of the 
Baltic Sea was shown in fig. 9. 
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Fig. 8. The isodynamic chart of the total magnetic 
flux pi of the Baltic Sea area 


Fig. 9. The real isodynamic chart of the total magnetic flux density 
of the Baltic Sea area (the same area as in fig.8) 


The variations of the magnetic field over the surface of the 
Baltic Sea are about 2 uT an area of several square kilometres. 
Such strong anomalies make serious problem in detection of 
ferromagnetic objects. 


MAGNETOMETERS 


The scalar magnetometers are used mainly in geological 
researches and military magnetic systems. The contemporary 
optically pumped magnetometers are of high sensitivity 107 T 
[8] (fig. 10). The SQUID vector magnetometers of the highest 
sensitivity (about 10°'* T) are also used [8, 9]. The magnetic 
field of the Earth is about 5:10% T (fig. 6). So, the sensitivity of 
the magnetometers is more than million fold less than the total 
magnetic field of the Earth. The magnetometers are installed 
on mobile platforms: planes, helicopters or UAV. 


MAGNETIC SYSTEMS ON AIR PLATFORMS 


The magnetic system installed on a plane has one ora set of 
magnetometers. The system with one magnetometer measures 
the total magnetic field (fig.11). 

On account of the natural magnetic variations (hundreds 
pT + several nT) the second magnetometer on the ground is 
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Fig. 10. The optically pumped magnetometer G-822A 
of Geometrics Company [13] 


Fig. 11. The military magnetic system with one magnetometer 
(the sensor is in the tail of the plane) 


used as the reference magnetometer. The system with a set 
of magnetometers measures the difference of modulus of the 
magnetic field [fig.12]. 


“4 E 


Fig. 12. The magnetic system with a set of magnetometers (the sensors 
are installed in the tail of the plane and on wings) [14] 


The magnetic system on the plane includes the digital 
compensator of magnetic disturbances [10, 11]. Each change 
ofthe position ofthe plane in respect of the vector of magnetic 
field of the Earth disturbs the magnetic field surrounding the 
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plane. These disturbances are caused by constant and induced 
magnetization of the ferromagnetic elements of the plane and 
eddy currents induced in the metallic elements of the plane. 
The mathematical model of the magnetic disturbances is used 
in the compensator [10]. The linear model is described by the 
equation: 


18 
H(t) =J k; AW (15) 
i=l 
where: 
k. — the constant factors of the model 


A(t) — the functions depend on the position of the plane 
t — time. 


The values of the functions A(t) are calculated on the 
basis of the direction cosines and the derivatives of the 
direction cosines [10]. The direction and derivative cosines 
are measured by 3-axis fluxgate magnetometer. The factors k, 
are connected with induced and constant magnetization of the 
ferromagnetic elements of the plane and with eddy currents. 
Every plane has some characteristic values of these factors. 
The digital compensator determines the factors k, on the basis 
of measurements of the magnetic flux density executed in the 
uniform magnetic field of the Earth with forced banking of the 
plane (pitch, yaw, roll). The simple magnetic system on a plane 
consists of a scalar magnetometer, a fluxgate magnetometer 
and a computer (fig. 13). 


fluxgate 
magnetometer 


scalar 
magnetometer 


Fig. 13. The structure of the magnetic system on a plane 


The example of the compensation of the magnetic 
disturbances of the plane was shown in fig.14. The magnetic 
system with a good compensator allows making measurements 
with sensitivity less than 100 pT. 
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Fig. 14. The results of the compensation of the disturbances 
of the magnetic field 


The first Polish navy aircraft with magnetic system was 
shown in fig.15. 

This patrol-pathfinder navy plane An—28B1R Bis is used 
for instance in searching and tracking of submarines in the 
Baltic Sea. The systems installed on a helicopter do not need the 


compensator, if the magnetic sensor is in the remote distance 
from the helicopter. The magnetic sensor is mainly inside the 
gondola, which is towed by the helicopter (fig.16). 


Fig. 15. The Polish navy aircraft with magnetic system 
(photo by: Mariusz Adamski) 


Fig. 16. The magnetic system on a helicopter 
with the magnetic sensor inside the gondola 


The example of the Polish navy helicopter Mi—14Pł with 


Fig. 17. The Polish navy helicopter Mi—14P! with magnetic system 
(photo by: Mariusz Adamski) 


DETECTION OF OBJECTS IN LOCAL 
ANOMALY AREAS 


The magnetic systems installed on planes or helicopters 
are used in geological researches and submarines or wracks 
detection. The strong local magnetic anomalies cause serious 
problems in detection of small signals connected with 
ferromagnetic objects. A good example of a special area with 


strong local anomalies is the Baltic Sea. The magnetic system 
should have a sophisticated detection algorithm. Particular 
problems occur in the system of a helicopter towing a gondola. 
The trajectory of the gondola’s flight is not a straight line. The 
deviations from the straight line of the gondola occur during 
the flight. The deviations are bigger during bad meteorological 
conditions (for example big crosswind). The gondola moves 
approximately along the trajectory: 


x(t) =v,tcos(a,)— Asin(27nf ,t)sin(a,) 


(16) 
y(t) =v,tsin(a,)+ Asin(2zf,t)cos(a,) 
where: 
v, — aconstant velocity of a platform 
A S amplitude of deviation 
f — frequency of deviation of the gondola 
a, — anangle between direction ofa straight line trajectory 


and x axis. 


The frequency of the gondola’s deviation depends on the 
length of the cable, mass and aerodynamics parameters of 
the gondola and meteorological conditions. The result of the 
gondola’s deviation is the interference signal in measurements 
of the magnetic flux density. The example of the local magnetic 
anomalies similar to real anomalies on the Baltic Sea was 
presented in fig. 18. 
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Fig. 18. The distribution of a local magnetic anomaly 
similar to real anomaly on the Baltic Sea 


The interference signal AB(t) appears in the measure signal 
with the frequency about f (fig.19): 


AB(t) = B x(t) - B yp (t) (17) 


where: 

B® — the changes of the magnetic flux density in time 
along a non-linear trajectory 

B y® — the changes of the magnetic flux density in time 
along a linear trajectory. 


The amplitude and waveband of the interference signal 
depend on many factors such as: 


e 


ko 


the distribution and values of the anomalies, 

the direction of the movements of the magnetic sensor, 
the mass of the gondola 

the length of the gondola’s cable 

the meteorological conditions. 


te 


ko 
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The values of these interferences could amount to several 
nT in extreme cases. 
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20 40 60 30 100 120 +140 160 180 200 
t [s] 
Fig. 19. The changes of the AB(t) caused by the gondola’ deviation 


(A=4m, a, = 0°, v= 50 m/s, f, = 0.1 Hz) 
along a linear trajectory a-b (fig.18) 


CONLUSIONS 


The measure methods, which are used in geological 
researches, are presented in this paper. The special attention 
was paid to magnetic systems installed on mobile platforms. 
The magnetic field of the Earth obtained from mathematical 
model was compared to the real magnetic field in the selected 
part of the Baltic Sea. The results of the calculations showed 
that the mathematical model of the Earth’s magnetic field 
does not consider local magnetic anomalies. The strong local 
magnetic anomalies cause serious problems with detection of 
underwater objects. Special problems appear in the magnetic 
system on a helicopter. The interference signals appear in the 
measurements of the magnetic flux density when the sensor is 
towed in the gondola. These signals have to be compensated 
in order to avoid false detection signals. 
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